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ABSTRACT 


Macrophages exist in most tissues and play a variety 
of functions in vertebrates. Teleost fish species are 
found in most aquatic environments throughout the 
world and are quite diverse for a group of vertebrate 
animals. Due to whole genome duplication and 
environmental adaptation, teleost — monocytes/ 
macrophages possess a variety of different functions 
and modulations compared with those of mammals. 
A deeper understanding of teleost | monocytes/ 
macrophages in the immune system will not only 
help develop teleost-specific methods of disease 
prevention but will also help improve our understanding 
of the various immune mechanisms in mammals. In 
this review, we summarize the differences in 
polarization and phagocytosis of teleost and 
mammalian macrophages to improve our 
understanding of the various immune mechanisms in 
vertebrates. 


Keywords: Teleost; | Monocytes /Macrophages; 
Phagocytosis; Cytokine production; Comparative 
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INTRODUCTION 


Macrophages exist in most animal tissues, in which they play 
a variety of functions. They are known by different names in 
various groups, such as amebocytes, hemocytes, coelomocytes, 
granulocytes, monocytes, and macrophages, but have similar 
morphology and comparable functions (Bilej et al., 2000; 
Wiegertjes et al., 2016). Macrophages are best known for 
their role in immunity, as elucidated by Eli Metchnikoff in the 
late nineteenth century (Tauber, 2003). Several studies have 
highlighted the wide range of functions of macrophages in 
vertebrate biology, including systemic metabolism, cold 
adaptation, tissue homeostasis, and development (Okabe & 
Medzhitov, 2016; Wynn et al., 2013). The basic functions of 
macrophages in vertebrates are cytokine production and 
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phagocytosis (Hodgkinson et al., 2015). In both mammals and 
fish, monocytes give rise to macrophages during inflammatory 
conditions (Hodgkinson et al., 2015). Moreover, macrophage 
colony stimulating factor (CSF1) plays a crucial role in 
macrophage growth and differentiation in both mammals and 
fish (Hodgkinson et al, 2015). Teleosts are among the 
evolutionarily oldest vertebrates, possessing both innate and 
classical adaptive immune systems (Dickerson & Findly, 
2017). Environmental factors have been shown to affect the 
immune system of teleosts (Makrinos & Bowden, 2016), which 
are widespread in most aquatic environments. Various immune 
genes in non-model teleosts have been identified, with 
transcriptome and genome development (Mackintosh & 
Ferrier, 2017; Qian et al., 2014; Shao et al., 2016), providing 
an opportunity to conduct studies on the adaptive evolution of 
the immune system. In this review, we focus on the differences 
in polarization and phagocytosis of teleost and mammalian 
macrophages, which should help in the development of a new 
perspective on macrophages and their role in adaptive evolution. 


DEFINITION OF 
TELEOSTS 


MONOCYTES/MACROPHAGES IN 


Monocytes/macrophages in mammals are an important 
component of the mononuclear phagocytic system, and play 
diverse roles during infection, inflammation, and tissue injury 
and repair (Okabe & Medzhitov, 2016). In mammals, protein 
markers are used to identify monocytes, macrophages, and 
dendritic cells by flow cytometry. Monocytes, which mainly 
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exist in bone marrow, blood, and spleen, can further differentiate 
into inflammatory macrophages and dendritic cells during 
inflammation (Ray & Rai, 2017; Shi & Pamer, 2011). Their 
migration to tissues and differentiation into inflammatory 
macrophages and dendritic cells are determined by the 
inflammatory milieu and pattern recognition receptors on the 
cells (Murray, 2018). Macrophages reside in a variety of 
tissues, including lymphoid and non-lymphoid tissues. They 
are equipped with a broad range of pattern recognition 
receptors to facilitate cytokine production and phagocytosis 
during inflammatory responses (Franz & Kagan, 2017). Classical 
dendritic cells are specialized antigen processing and presenting 
cells that exhibit high phagocytic activity as immature cells 
and high cytokine-producing capacity as mature cells (Qian & 
Cao, 2018). 

In teleosts, monocytes/macrophages have been identified in 
a variety of species, including rainbow trout (Oncorhynchus 
mykiss) (Leal et al., 2017), zebrafish (Danio rerio) (Yu et al., 
2017), goldfish (Carassius auratus L.) (Hodgkinson et al., 
2017), and ayu (Plecoglossus altivelis) (Lu et al., 2017). The 
human CD14 antigen is highly expressed in monocytes and to 
a lesser extent in macrophages (Shi & Pamer, 2011; Zhang et 
al., 2017). However, CD14 has not been identified in the 
teleost genome (Novoa et al., 2009). Hence, it is difficult to 
discriminate between teleost monocytes and macrophages 
based on protein markers using flow cytometry. Stafford et al. 
(2001) developed a method using cell-conditioned medium to 
obtain monocyte-like cells in rainbow trout. However, it is 
unclear whether monocytes can be identified from the culture 
of teleost monocytes/macrophages based on acute isolation 
using simple procedures, such as flow cytometry. Although 
dendritic-like cells have been identified in rainbow trout (Soleto 
et al., 2018), it is unclear whether such cells exist in other 
dominant teleost species. As monocytes and macrophages 
share similar characteristics, most teleost investigations have 
defined the adhered mononuclear cells of the kidney or head 
kidney as monocytes/macrophages or macrophages. 

Phagocytosis is the defining characteristic of macrophages 
(large eaters), as classified by Metchnikoff (Tauber, 2003). 
However, most phagocytes among the peripheral blood 
mononuclear cells of teleosts are IgM* B lymphocytes, not 
monocytes/macrophages (Li et al., 2006). Thus, morphological 
analysis, apart from phagocytosis index analysis, is necessary 
to identify the monocytes/macrophages in teleosts. 


MACROPHAGE 
PRODUCTION 


POLARIZATION AND CYTOKINE 


In vertebrates, the inflammatory response of macrophages 
plays an important role against pathogens (Geissmann et al., 
2010; Ricci et al., 2018). Macrophage polarization against 
different inflammatory stimuli depends on environmental cues 
or pathophysiological conditions (Lawrence & Natoli, 2011). 
Classically activated macrophages (M1 type) are induced by 
lipopolysaccharides (LPSs), a major component of the outer 
membrane of gram-negative bacteria, and IFN-y, and express 


pro-inflammatory mediators. Conversely, alternatively activated 
macrophages (M2 type) are induced by IL-4 and IL-13, and 
express high levels of anti-inflammatory mediators (Shapouri- 
Moghaddam et al., 2018). Macrophage polarization is also 
regulated by soluble proteins, intracellular signals, and 
transcription factors. Galectin-dependent regulatory signaling 
stimulates M2-type macrophage polarization (Blidner et al., 
2015). Toll-like receptor (TLR) signaling activates the signal 
transducer and activator of transcription (STAT) 1 protein to 
skew macrophage function towards the M1 type, whereas 
activation of STAT3 by IL-4 and IL-13 skews macrophage 
function towards the M2 type (Sica & Mantovani, 2012). 
Similarly, the ablation of protein kinase B a (PKBa/Akt1) or 
protein kinase B B (PKBf/AKkt2) differentially affect macrophage 
polarization (Arranz et al., 2012). Tissue milieus with M1 type 
are biased towards cell-mediated cytotoxicity, whereas the 
term "M2 type" is used for a variety of conditions that inhibit 
M1 type (Yamaguchi et al., 2015). The immune milieus are 
skewed to M2 type in some tissues, like the gills in teleosts 
and uterus in pregnant mammals (Yamaguchi et al., 2015). 
Monocyte/macrophage polarization has also been detected 
in teleosts (Wiegertjes et al., 2016). Cytokines participate in 
teleost monocyte/macrophage polarization, particularly IFN - y 
and IL-4/IL-13 (Wiegertjes et al., 2016). Moreover, inducible 
nitric oxide synthase (iNOS) is a marker for M1 type and 
arginase 2 is a marker for M2 type in teleost monocytes / 
macrophages (Wiegertjes et al., 2016). In carp, the expression 
of pro-inflammatory cytokines IL-18, TNF-a, and CXCa, peak 
in peritoneal leukocytes 24 h after zymosan induction, 
whereas the expression of anti-inflammatory mediators IL-10 
and arginase 2 peak 96 h and 168 h after zymosan induction, 
respectively (Chadzinska et al., 2008). This suggests that 
monocytes/macrophages display both classic- and alternative 
pathway-induced polarization upon immune stimulation in 
vivo. Мопосуѓе / macrophage polarization in teleosts has also 
been investigated in vitro. LPSs from gram-negative bacteria 
are probably the best studied microbial stimuli for macrophage 
activation. However, the mammalian LPS receptor, TLR4, may 
not be functional in teleosts (Sepulcre et al., 2009). The 
presence of LPSs may be sensed by other mechanisms in 
teleosts, as LPSs are still thought to be important immune 
stimulators (Meng et al., 2012). Nitrite production in in vitro 
monocyte/macrophage culture is up-regulated after LPS 
stimulation in teleosts (Joerink et al., 2006). Other pro- 
inflammatory mediators, such as IL-18, IL-12, and iNOS are 
also up-regulated in LPS-stimulated teleost monocytes/ 
macrophages (Joerink et al., 2006), suggesting that LPS 
induces an М1-уре polarization in teleost monocytes/ 
macrophages. In mammals, M2-type polarization is mainly 
induced by IL-4, IL-13, parasite infection, CSF1, TGF-B, and 
IL-10 (Sica & Mantovani, 2012). Apart from anti-inflammatory 
mediators, cAMP plays an important role in M2-type monocyte/ 
macrophage polarization signaling (Bystrom et al., 2008). In 
teleosts, at least two IL-4/IL-13 genes exist (IL-4/13A and IL-4/ 
13B), both with low homology to IL-4 and IL-13 (Ohtani et al., 
2008). In goldfish, recombinant IL-4/13 has been found to 
induce arginase activity and down-regulate the nitric oxide 
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(NO) response in primary 
(Hodgkinson et al., 2017), suggesting IL-4/13 functions to 
induce | M2-type polarization in teleost monocytes / 
macrophages. In teleost carp (Joerink et al., 2006) and ayu 
(Chen et al., 2018), cAMP has also been employed to 
successfully induce M2-type monocyte/macrophage 
polarization. 

Several new mechanisms for macrophage polarization have 
been identified recently. Mammals possess single IFN-y 
molecules that are important in the activation of M1-type 
macrophages (Grayfer et al., 2018). Many teleosts have 
multiple IFN-y molecules, some of which are elicitors of 
reactive oxygen species (ROS) but not NO, whereas others 
elicit NO but limited ROS (Grayfer et al., 2018). Although 
functional analogs of the mammalian M1/M2 macrophage 
subsets are present in teleosts, defining the regulatory 
mechanisms governing the polarization of these effector 
populations is a far more challenging goal (Grayfer et al., 
2018). Many immune genes exist in two copies in the teleost 
genome due to genome duplication (Aghaallaei et al., 2010). 
These redundant genes may result in sub-functionalization, as 
in the case of European sea bass (Dicentrarchus labrax), in 
which different hepcidins exhibit different roles (Neves et al., 
2015). We observed that ayu has two CXCR3 genes, which 
contribute to monocyte / macrophage polarization (Lu et al., 
2017). In mammals, the chemokine receptor CXCR3 exists as 
a single gene, and is preferentially expressed on immune cells 
to aid in cell migration to the sites of inflammation (Bromley et 
al., 2008). Furthermore, CXCR3.1 (CXCR3b) and CXCR3.2 
(CXCR3a) are found in zebrafish (Danio rerio), Japanese 
ricefish (Oryzias latipes), spotted green pufferfish (Tetraodon 
nigroviridis), ayu (Plecoglossus altivelis), and grass carp 
(Ctenopharyngodon idella) (Aghaallaei et al., 2010; Lu et al., 
2017). However, more research is necessary to illustrate the 
teleost-specific mechanisms of monocyte/ macrophage 
polarization. 


monocytes/macrophages 


PHAGOCYTOSIS BY MONOCYTES/MACROPHAGES 


Phagocytosis is an important cellular process for the induction 
of antimicrobial responses and regulation of adaptive immunity 
(Rieger et al., 2012). After phagocytosis, both teleost and 
mammalian macrophages show pro-inflammatory апа 
homeostatic responses (Rieger et al., 2012). Mammalian 
neutrophils have the capacity to internalize apoptotic bodies, 
whereas teleost neutrophils do not possess the same activity 
(Rieger et al., 2012). Most studies have shown that in fish, 
monocytes, macrophages, and neutrophils are the main 
phagocytic cells, as found in mammals (Esteban et al., 2015). 
Furthermore, phagocytosis in teleosts has been observed in 
several other kinds of cells, including B-1 cells (Li et al., 
2006),  yó-T cells (Wan et al., 2016), and thrombocytes 
(Nagasawa et al., 2014). Although mammalian phagocytosis 
has also been observed in B-1 cells (Parra et al., 2012) and 
y5-T cells (Wu et al., 2009), its presence in thrombocytes 
seems unique to teleosts. It is well established that teleost 
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monocytes/macrophages generate reactive oxygen 
intermediates as part of antimicrobial defense, similar to that 
observed in mammalian macrophages (Hodgkinson et al., 
2015). However, significant differences between goldfish and 
mice have been observed in response to pro-inflammatory 
and homeostatic internalization signals during phagocytosis of 
immune cells (Rieger et al., 2012). When European sea bass 
(Dicentrarchus labrax) monocytes/macrophages have been 
incubated with different pathogenic agents, different pathogens 
have been observed to have different effects on monocyte/ 
macrophage activity (Bennani et al., 1995). However, further 
investigation is necessary to illustrate the specific mechanisms 
of pathogen recognition and phagocytosis in teleosts. 

Phagocytosis in mammals is triggered by pathogen 
recognition, which is a complex process involving a variety of 
pattern-recognition receptors. The main pattern-recognition 
receptors include lectin-like recognition molecules, C-type 
lectins, scavenger receptors, mannose receptors, CD14, and 
Toll-like receptors (Uribe-Querol & Rosales, 2017). In teleosts, 
the pathogen recognition mechanism is different from that of 
mammals. LPSs are sensed by a variety of pattern-recognition 
receptors (Ranf, 2016) present on the cell surface (Triantafilou 
& Triantafilou, 2002). Phagocytosis in macrophages is regulated 
by LPS recognition receptors, such as TLR4 (Lv et al., 2017) 
and CD14 (Lingnau et al., 2007). TLR4 was the first receptor 
identified to be involved in the recognition of LPSs. TLR4 must 
be associated with myeloid differentiation protein 2 (MD2) for 
interaction with LPSs, and activation of TLR4 is preceded by 
the transfer of LPSs to CD14 by an LPS-binding protein 
(Neyen & Lemaitre, 2016). However, CD14 and MD2 do not 
exist in teleost genomes (lliev et al., 2005). Thus, it remains 
unclear how teleost monocytes/macrophages compensate for 
the CD14 function during phagocytosis of pathogens. 


FUTURE DEVELOPMENTS 


Teleost fish species are found throughout the world and are 
quite diverse for a group of vertebrate animals. Monocytes/ 
macrophages are the basic immune cells not only in mammals 
but also in teleosts. Monocytes/macrophages are easy to 
isolate, purify, and segregate, providing us with a useful tool 
for understanding the differences between the immune 
Systems of mammals and teleosts. Although the basic function 
of monocytes/macrophages is similar in vertebrates, there are 
a variety of different points between mammals and teleosts 
(Figure 1). These differences have likely arisen from genetic 
factors, such as whole genome duplication, and environmental 
adaptations. Whole genome duplication in teleosts produces a 
variety of redundant genes, which may be sub-functional 
during adaptive evolution. Gene function in teleosts is 
modulated by environmental factors, such as salinity, hypoxic 
conditions, and temperature. Therefore, future investigations 
regarding the mechanisms of teleost monocyte / macrophage 
polarization and phagocytosis under the influence of various 
environmental factors are necessary. A deeper understanding 
of the teleost immune system will not only help develop 


teleost-specific methods of disease prevention but will also 
help improve our understanding of the various immune 
mechanisms within mammals. 


Teleosts 


| 


Mammals 
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ABSTRACT 


We describe a new species of the genus Tylototriton 
from Ingyin Taung Mt., Mohnyin Township, Kachin 
State, Myanmar, based on morphological and 
molecular evidence. The new species is assigned to 
the subgenus Tylototriton s. str. and is clearly distinct 
from all known congeners by the following characters: 
medium body size; thin, long tail, lacking lateral 
grooves; rough skin; truncate snout; wide, protruding 
supratemporal bony ridges on head, beginning at 
anterior corner of orbit; weak, almost indistinct sagittal 
ridge; long, thin limbs, broadly overlapping when 
adpressed along body; distinct, wide, non-segmented 
vertebral ridge; 13 or 14 rib nodules; brown to dark- 
brown background coloration with dull orange-brown 
to yellowish-brown markings оп labial regions, 
parotoids, rib nodules, whole limbs, vent, and ventral 
tail ridge. We also briefly discuss biogeography and 
species diversity of the genus Tylototriton in Myanmar. 


Keywords: Tylototriton kachinorum sp. nov.; mtDNA 
genealogy; ND2; 16S rRNA; Shan; Biogeography; 
Endemism; Taxonomy 


INTRODUCTION 


The salamandrid genus Tylototriton Anderson, 1871, or 
Crocodile Newts, currently includes 24 recognized species, 
inhabiting montane forest areas throughout the Asian 
monsoon climate zone from eastern Himalaya, southern and 
central China including Hainan Island, to northern Indochina 
including Vietnam, Laos, Thailand, and Myanmar (Hernandez, 
2016; Wang et al., 2018). Recent progress in phylogenetic 
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studies of the genus Tylototriton has indicated that the genus 
is monophyletic (Nishikawa et al., 2013a, 2013b; Phimmachak 
et al., 2015; Wang et al., 2018) and includes two major 
groups, corresponding to the subgenera Tylototriton s. str. and 
Yaotriton (Wang et al., 2018). Molecular taxonomy methods 
have proven to be useful for deciphering taxonomic diversity 
of the genus Tylototriton, with 13 species (over 50%) described 
in the past five years (Grismer et al., 2018a; Hou et al., 2012; 
Khatiwada et al., 2015; Nishikawa et al., 2013a, 2013b, 2014; 
Phimmachak et al., 2015; Shen et al., 2012; Yang et al., 2014; 
Zhao et al., 2012). However, recent molecular surveys indicate 
that our knowledge on taxonomic diversity of the genus 
Tylototriton is still far from complete, revealing several cryptic 
lineages likely corresponding to as yet undescribed species 
(Grismer et al., 2018a; Wang et al., 2018). 

Myanmar, previously known as Burma, is the largest country 
of mainland Southeast Asia. Despite this, its herpetofauna 
remains one of the least explored in the region (Grismer et al., 
2018a). Members of the genus Tylototriton have long been 
recorded from northern and eastern parts of Myanmar and 
have been traditionally classified as T. verrucosus Anderson, 
1871 (Gyi, 1969). Nishikawa et al. (2014), based on the 
examination of specimens assigned to T. verrucosus collected 
from the Shan Plateau in eastern Myanmar and pet-trade 
animals assumed to originate from Myanmar, recently described 
a new species, 7T. shanorum Nishikawa, Matsui & Rao, 2014. 
Soon after, Phimmachak et al. (2015) published sequence 
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data for specimens collected from the Sagaing Region and 
Kachin State in northern Myanmar, which were reported as T. 
verrucosus. More recently, Grismer et al. (2018a) demonstrated 
the presence of two morphologically and genetically distinct 
lineages of Tylototriton in the Shan Plateau and described a 
new species from its north-western edge, T. ngarsuensis 
Grismer, Wood, Quah, Thura, Espinoza, Grismer, Murdoch & 
Lin, 2018. The same work of Grismer et al. (2018a) reanalyzed 
sequences of specimens from the Sagaing and Kachin 
regions of Myanmar and demonstrated that they belong to a 
distinct lineage — Tylototriton sp. 1, distinct from T. verrucosus 
s. str. Recent work also indicated the presence of T. 
himalayanus Khatiwada, Wang, Ghimire, Vasudevan, Paudel 
& Jiang, 2015 (a species described from Nepalese Himalaya) 
in northern Myanmar, but without providing voucher specimen 
information or any other justification for this identification 
(Hernandez, 2016; Hernandez et al., 2018). The recent 
monographic review of the genus Tylototriton by Hernandez 
(2016) also indicated the possibility of the occurrence of T. 
uyenoi Nishikawa, Khonsue, Pomchote & Matsui, 2013 and T. 
shanjing Nussbaum, Brodie & Yang, 1995 in parts of Myanmar 
adjacent to northern Thailand and the southwestern Yunnan 
Province of China; however, these records are not supported 
by voucher specimens. Thus, our knowledge on the taxonomic 
composition and diversity of the genus Tylototriton in Myanmar 
is still far from complete. 


1 
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In the present paper, we report on a new population of the 
genus Tylototriton from the Kachin Hills in the southern part of 
Kachin State, northern Myanmar. We applied morphological 
and molecular methods to evaluate its taxonomic status and 
describe it as a new species. We also discuss biogeography 
and taxonomy of the genus Tylototriton in Myanmar. 


MATERIALS AND METHODS 


Sample collection 

Fieldwork was carried out in northern Myanmar, Kachin State, 
from 14 to 21 July 2018. Specimens of Tylototriton sp. were 
collected by hand in swamps in forest clearings surrounded by 
montane evergreen tropical forests of Ingyin Taung Mountain, 
Indawgyi Lake area, Kachin State (Figure 1; samples 22-23). 
Geographic coordinates and altitude were obtained using a 
Garmin GPSMAP 60CSx GPS receiver (Garmin Ltd., USA) 
and recorded in datum WGS 84. Specimens were euthanized 
by 2096 benzocaine and tissue samples for genetic analysis 
were taken and stored in 9696 ethanol (femoral muscles) prior 
to preservation. Specimens were subsequently preserved in 
7096 ethanol and deposited in the herpetological collection of 
the Zoological Museum of Moscow State University (ZMMU) 
in Moscow, Russia; Zoological Institute of Russian Academy 
of Sciences (ZISP) in St. Petersburg, Russia; and Zoology 
Department of University of Mandalay (ZDUM), Mandalay, 
Myanmar. 


Figure 1 Distribution of Tylototriton and sampling localities examined in this study 


Colors of dots correspond to clades 1—5 denoted in Figure 2. Sample numbers 1-115 


kachinorum sp. nov. is marked with a star. Photo by Nikolay A. Poyarkov. 
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refer to Table 1. Type locality of new species Tylototriton 


Morphological description 

Specimens of Tylototriton sp. were photographed in life and 
after preservation. The sex and maturity of the specimens and 
number of eggs were checked and counted by minor 
dissections. Measurements were taken using a digital caliper 
to the nearest 0.01 mm, subsequently rounded to 0.1 mm. We 
used a_ stereoscopic light binocular microscope when 
necessary. Statistical analyses were performed with Statistica 
8.0 (Version 8.0; StatSoft, Tulsa, OK, USA). 


Adult morphology 
Morphometrics followed Nishikawa et al. (2014), Khatiwada et 
al. (2015), and Okamiya et al. (2018) and included the following 
24 measurements: (1) SVL (snout-vent length) from tip of 
snout to anterior tip of vent; (2) HL (head length); (3) HW 
(head width); (4) MXHW (maximum head width); (5) IND 
(internarial distance); (6) AGD (axilla-groin distance); (7) TRL 
(trunk length); (8) TAL (tail length) from anterior tip of vent to 
tail tip; (9) VL (vent length); (10) FLL (forelimb length); (11) HLL 
(hindlimb length); (12) VTW (vomerine tooth series width): 
greatest width of vomerine tooth series; (13) VTL (vomerine 
tooth series length): greatest length of vomerine tooth series; 
(14) LJL (lower jaw length from tip of lower jaw to articulation 
of upper and lower jaws); (15) SL (snout length from tip of 
snout to anterior tip of upper eyelid); (16) IOD (minimum 
interorbital distance); (17) UEW (maximum upper eyelid width); 
(18) UEL (upper eyelid length, distance between anterior and 
posterior angles); (19) OL (orbit length); (20) BTAW (basal tail 
width at level of anterior tip of cloaca); (21) MTAW (tail width at 
mid-level of tail); (22) MXTAH (maximum tail height); (23) MTAH 
(tail height at mid-level of tail); (24) ON (orbitonarial distance). 
For holotype description, we examined the following 12 
morphometric and three meristic characters following Poyarkov 
et al. (2012) and Okamiya et al. (2018): additional morphometric 
characters: (25) ICD (intercanthal distance); (26) CW (chest 
width); (27) NSD (nostril-snout distance); (28) 1FL (first finger 
length from base to tip); (29) 2FL (second finger length from base 
to tip); (30) 3FL (third finger length from base to tip); (31) 4FL 
(fourth finger length from base to tip); (32) 1TL (first toe length 
from base to tip); (33) 2TL (second toe length from base to 
tip); (34) 3TL (third toe length from base to tip); (35) 4TL (fourth 
toe length from base to tip); (36) 5TL (fifth toe length from 
base to tip); meristic characters: (37) UJTN (number of teeth 
onupper jaw); (38) LUTN (number ofteeth onlower jaw); (39) VTN, 
number of teeth on vomer. We also recorded the following 
characters as per Nishikawa et al. (2014) and Grismer et al. 
(2018a): shape of vomerine teeth and their positional 
relationship relative to choanae; skin texture; number and 
shape of rib nodules counted from posterior margin of vent to 
axilla; width and prominence of vertebral ridge and head 
ridges; and coloration of dorsum, venter, head, labial region, 
parotoid glands, rib nodules, limbs, soles, palms, tail surfaces, 
and vent region. 

The characters of adult morphology chosen for 
comparison and data on other Tylototriton species were 
taken from the following sources: Anderson (1871); Bóhme 


et al. (2005); Chen et al. (2010); Fang & Chang (1932); Fei 
et al. (1984); Grismer et al. (2018a); Hernandez (2016); Hou 
et al. (2012); Khatiwada et al. (2015); Le et al. (2015); Liu 
(1950); Nishikawa et al. (2013a; 2013b; 2014); Nussbaum et 
al. (1995); Phimmachak et al. (2015); Qian et al. (2017); 
Shen et al. (2012); Stuart et al. (2010); Unterstein (1930); 
Yang et al. (2014). 


Larval morphology 

Description of larval morphology followed Okamiya et al. 
(2018). For larval specimens, we recorded nine morphometric 
characters including (1) SVL, (2) HL, (3) HW, (4) OL, (5) AGD, 
(6) TAL, (7) FLL, (8) HLL, and (9) MXTAH (definition same as 
for adult morphology). Developmental stages were determined 
following Grosse (2013). 


DNA isolation, PCR, and sequencing 

Total genomic DNA was extracted from 9596 ethanol-preserved 
muscle tissues using standard phenol-chloroform extraction 
protocols (Hillis et al., 1996). Total DNA concentration was 
estimated in 1 uL using a NanoDrop 2000 spectrophotometer 
(Thermo Scientific, USA), and consequently adjusted to 
100 ng DNA/LL. 

We amplified two mtDNA fragments consisting of partial 
sequences of the ND2 and 16S rRNA mtDNA genes. These 
markers were chosen as they are useful in studies of Tylototriton 
phylogeny and taxonomy (Nishikawa et al., 2013a, 2013b, 
2014; Wang et al., 2018; Zhao et al., 2012 and references 
therein). We used the 16L-1 (forward) (5-CTGACCGTGCAAA 
GGTAGCGTAATCACT-3') and 16H-1 (reverse) (5'- CTCCGG 
TCTGAACTCAGATCACGTAGG-3') primers to amplify the 
16S rRNA fragments (Hedges, 1994). For amplification and 
sequencing of the ND2 gene, we used the SL-1 (forward) (5'- 
ATAGAGGTTCAAACCCTCTC-3) and SL-2 (reverse) (5'- 
TTAAAGTGTCTGGGTTGCATTCAG-3') primers of Wang et 
al. (2018). Polymerase chain reaction (PCR) was performed 
in 20 pL reactions using 50 ng genomic DNA, 10 nmol of 
each primer, 15 nmol of each dNTP, 50 nmol additional 
MgCL, Taq PCR buffer (10 mmol/L Tris-HCl, pH 8.3, 
50 mmol/L KCI, 1.1 mmol/L М9СІ,, and 0.01% gelatin), and 
1 U of Taq DNA polymerase. PCR cycles included an initial 
denaturation step of 4 min at 94 °C and 35 cycles of 
denaturation for 30 s at 94 °C, primer annealing for 30 s at 
48—58 °C, and extension for 1 min 30 s at 72 °C. PCR 
products were visualized by agarose gel electrophoresis in 
the presence of ethidium bromide and consequently purified 
using 2 uL from a 1:4 dilution of ExoSaplt (Amersham, UK) 
per 5 uL of PCR product prior to cycle sequencing. 
Sequencing was performed in both directions using the same 
primers as used in PCR on an ABI3730xl automated sequencer 
(Applied Biosystems, USA) at Evrogen Inc., Moscow (Russia). 
The newly obtained sequences were aligned and deposited 
in GenBank under the accession numbers MK095616— 
MK095619 and MK097271—MK097274 (Table 1). Sequences 
of 25 other Tylototriton species used for comparisons were 
obtained from GenBank (Table 1). 
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Table 1 


Sample 


Sequences and voucher specimens of Tylototriton and outgroup taxa used in this study 





No. Species name Voucher number Locality 16S rRNA ND2 
Ingroup: 

1 Tylototriton taliangensis CIB GG200110183 Shimian Co., Ya'an City, Sichuan, China KY800559  KC147819 
2 Tylototriton taliangensis CIB GG200110185 Shimian Co., Ya'an City, Sichuan, China KY800560  KY800829 
3 Tylototriton taliangensis CIB GG200110186 Shimian Co. Ya'an City, Sichuan, China KY800561 КҮ800830 
4 Tylototriton pseudoverrucosus СІВ WCG2012003 ИА me Yi Autonomous Prefec- к увоо5от KY800861 
5 Tylototriton pseudoverrucosus СІВ WCG2012007 AAR a е Yi Autonomous Prefec- с увоо5ов күв00862 
6 Tylototriton pseudoverrucosus CIB WCG2012012 A ma Yi Autonomous Prefec" Kvgogs99  KY800860 
7 Tylototriton kweichowensis СІВ WG20080818014 Віјіе City, Guizhou, China KY800551 KY800823 
8 Tylototriton kweichowensis СІВ WG20080818018 Віјіе City, Guizhou, China KY800552 КҮ800824 
9 Tylototriton kweichowensis CIB 20050213 Shuicheng Co., Guizhou, China KY800557  KY800827 
10 Tylototriton kweichowensis CIB 20050215 Shuicheng Co., Guizhou, China KY800558  KY800828 
11 Tylototriton shanorum lineage 1 CAS 230933 Taunggyi Township, Shan, Myanmar = AB922822 
12 Tylototriton shanorum lineage 1 CAS 230940 Taunggyi Township, Shan, Myanmar = AB922823 
13 Tylototriton shanorum lineage 1 CAS 230899 Taunggyi Township, Shan, Myanmar = HM770087 
14 Tylototriton shanorum lineage 2 KUHE 42348 pet trade, presumably from Myanmar = AB769544 
15 Tylototriton ngarsuensis LSUHC 13762 Ngar Su, Shan, Myanmar = MH836585 
16 Tylototriton ngarsuensis LSUHC 13763 Ngar Su, Shan, Myanmar = MH836584 
17 Tylototriton sp. 1 CAS 245418 Chipwi Township, Kachin State, Myanmar = KT304279 
18 Tylototriton sp. 1 CAS 245290 Lahe Township, Sagaing Region, Myanmar == KT304278 
19 Tylototriton himalayanus CIB 201406246 Mai Pokhari, Illam, Mechi, Nepal KY800590 KT765173 
20 Tylototriton himalayanus CIB 201406284 Mai Pokhari, Illam, Mechi, Nepal KY800591 KT765207 
21 Tylototriton himalayanus CIB 201406285 Mai Pokhari, Illam, Mechi, Nepal KY800592 KT765208 
22 Tylototriton kachinorum sp. nov. ZMMU А5953 Ingyin Taung Mt., Indawgyi, Kachin, Myanmar MK095618 МК097273 
23 Tylototriton kachinorum sp. nov. ZDUM-0103 Ingyin Taung Mt., Indawgyi, Kachin, Myanmar MK095619 — MK097274 
24 Tylototriton yangi KUHE 42282 Pingbian Co., Yunnan, China KY800624 KY800887 
25 Tylototriton uyenoi ZMMU NAP-08220 Doi Inthanon, Chiang Mai, Thailand MK095617 — MK097272 
26 Tylototriton uyenoi KUHE 19037 Doi Inthanon, Chiang Mai, Thailand = AB830730 
27 Tylototriton uyenoi KUHE 19147 Doi Ang Khang, Chiang Mai, Thailand — AB830729 
28 Tylototriton uyenoi ITNT1 Doi Suthep, Chiang Mai, Thailand = AB830733 
29 Tylototriton anguliceps TBU PAE671 Thuan Chau, Son La, Vietnam — LC017833 
30 Tylototriton anguliceps VNMN A20143 Muong Nhe, Dien Bien, Vietnam == LC017832 
31 Tylototriton anguliceps LK003 Doi Lahnga, Chiang Mai, Thailand — AB830728 
32 Tylototriton podichthys KUHE 34399 Phu Pan, Xam Neua, Laos == АВ830727 
33 Tylototriton podichthys IEBR A2014-1 Xam Neua, Huaphanh, Laos = LC017835 
34 Tylototriton pulcherrimus CIB TY 040 Lüchun Co., Yunnan, China KY800626 KY800890 
35 Tylototriton pulcherrimus KUHE 46406 Pet Trade KY800620 КҮ800880 
36 Tylototriton verrucosus CIB TSHS1 Longchuan Co., Dehong, Yunnan, China KY800581 KY800847 
ЭТ Tylototriton verrucosus CIB-TSHS2 Longchuan Co., Dehong, Yunnan, China KY800582 KY800848 
38 Tylototriton verrucosus CIB TSHS3 Longchuan Co., Dehong, Yunnan, China KY800583 KY800849 
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Sample 


Continued 





No. Species name Voucher number Locality 16S rRNA ND2 

39 Tylototriton verrucosus CIB TSHS4 Longchuan Co., Dehong, Yunnan, China KY800584 KY800850 
40 Tylototriton verrucosus CIB TSHS5 Longchuan Co., Dehong, Yunnan, China KY800585 КҮ800851 
41 Tylototriton verrucosus CIB TSHS6 Longchuan Co., Dehong, Yunnan, China KY800586  KY800852 
42 Tylototriton shanjing KIZ 201306081 Yongde Co., Yunnan, China KY800593 КҮ800856 
43 Tylototriton shanjing KIZ 201306098 Yun Co., Yunnan, China KY800594 КҮ800857 
44 Tylototriton shanjing KIZ 201306102 Jingdong Co., Yunnan, China KY800595 КҮ800858 
45 Tylototriton shanjing KIZ 201306108 Nanjian Co., Yunnan, China KY800596 КҮ800859 
46 Tylototriton shanjing CIB 980004 Baoshan City, Yunnan, China KY800562 КҮ800831 
47 Tylototriton shanjing CIB 980005 Baoshan City, Yunnan, China KY800563 КҮ800832 
48 Tylototriton shanjing CIB 980006 Baoshan City, Yunnan, China KY800564 КҮ800833 
49 Tylototriton panhai ZMMU NAP-08217 Phu Hin Rong Kla NP, Phitsanulok, Thailand MK095616 | MK097271 
50 Tylototriton panhai PH019 Phu Hin Rong Kla NP, Phitsanulok, Thailand = AB830735 
51 Tylototriton panhai PLOO9 Phu Luang WS, Loei, Thailand = AB830736 
52 Tylototriton vietnamensis IEBR A.3674 Yen Tu, Bac Giang, Vietnam KY800614 KY800874 
53 Tylototriton vietnamensis IEBR A.0702 Mauson, Lang Son, Vietnam KY800612 KY800872 
54 Tylototriton vietnamensis IEBR A.0701 Mauson, Lang Son, Vietnam KY800613 KY800873 
55 Tylototriton ziegleri VNMN 3389 Bao Lac, Cao Bang, Vietnam = KY800888 
56 Tylototriton ziegleri VNMN 3390 Quan Ba, Ha Giang, Vietnam KY800625 КҮ800889 
57 Tylototriton ziegleri VNUH HG.082 Quan Ba, Ha Giang, Vietnam KY800610 KY800870 
58 Tylototriton ziegleri VNUH HG.081 Quan Ba, Ha Giang, Vietnam KY800611 КҮ800871 
59 Tylototriton hainanensis CIB 20081048 Mt. Diaoluo, Hainan, China KY800553 KC147817 
60 Tylototriton hainanensis CIB 20081049 Mt. Diaoluo, Hainan, China KY800554 KC147818 
61 Tylototriton hainanensis CIB 20081051 Mt. Diaoluo, Hainan, China KY800555  KY800825 
62 Tylototriton hainanensis CIB 20081052 Mt. Diaoluo, Hainan, China KY800556 КҮ800826 
63 Tylototriton asperrimus lineage 2 CIB XZ20091201 Xinyi City, Guangdong, China KY800616 КҮ800876 
64 Tylototriton asperrimus lineage 2 CIB XZ20091204 Xinyi City, Guangdong, China KY800619 КҮ800879 
65 Tylototriton asperrimus lineage 2 CIB XZ20091 Xinyi City, Guangdong, China KY800618 КҮ800878 
66 Tylototriton asperrimus lineage 2 CIB XZ20092 Xinyi City, Guangdong, China KY800617 КҮ800877 
67 Tylototriton notialis FMNH HERP271120 Boualapha Dist., Khammouan, Laos = HM462061 
68 Tylototriton notialis FMNH HERP271121 Boualapha Dist., Khammouan, Laos — HM462062 
69 Tylototriton notialis VNMN TAO1229 Pu Hoat, Nghe An, Vietnam — KY800883 
70 Tylototriton notialis VNMN TAO1235 Pu Hoat, Nghe An, Vietnam = KY800884 
71 Tylototriton asperrimus lineage 1 ММММ TAO1213 Thuong Tien, Hoa Binh, Vietnam KY800623 КҮ800885 
72 Tylototriton asperrimus lineage 1 ММММ TAO1214 Thuong Tien, Hoa Binh, Vietnam = KY800886 
73 Tylototriton asperrimus lineage 1 СІВ 70063 Longsheng Co., Guangxi, China KY800549 КС147816 
74 Tylototriton asperrimus lineage 1 СІВ GX20080714 Jinxiu Co., Guangxi, China KY800546 КҮ800819 
75 Tylototriton asperrimus lineage 1 CIB GX200807010 Jinxiu Co., Guangxi, China KY800547 КҮ800820 
76 Tylototriton asperrimus lineage 1 СІВ GX200807012 Jinxiu Co., Guangxi, China KY800548 KY800821 
77 Tylototriton asperrimus lineage 1 СІВ GX200807016 Jinxiu Co., Guangxi, China KY800550 КҮ800822 
78 Tylototriton asperrimus lineage 1 CIB 20070715 Jinxiu Co., Guangxi, China KY800565 КҮ800834 
79 Tylototriton asperrimus lineage 1 CIB 200807055 Jinxiu Co., Guangxi, China KY800566 КС147815 
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Continued 








No Species name Voucher number Locality 16S rRNA ND2 

80 Tylototriton liuyangensis CSUFT 20100108 Liuyang City, Hunan, China KY800606 КЈ205598 
81 Tylototriton liuyangensis CIB 110601F06 Liuyang City, Hunan, China KY800615 КҮ800875 
82 Tylototriton lizhenchangi KUHE 42316 Yizhang Co., Hunan, China KY800621 KY800881 
83 Tylototriton lizhenchangi KUHE 42317 Yizhang Co., Hunan, China KY800622 KY800882 
84 Tylototriton dabienicus lineage 2 СІВ 08042905-2 Yuexi Co. Anhui, China KY800587 КҮ800853 
85 Tylototriton dabienicus lineage 2 СІВ 08042905-3 Yuexi Co. Anhui, China KY800588 КҮ800854 
86 Tylototriton dabienicus lineage 2 СІВ 08042905-4 Yuexi Co. Anhui, China KY800589 КҮ800855 
87 Tylototriton broadoridgus CIB 200085 Sangzhi Co., Hunan, China KY800569 КС147814 
88 Tylototriton broadoridgus CIB 200084 Sangzhi Co., Hunan, China KY800570  KY800837 
89 Tylototriton dabienicus lineage 1 HNNU 1004-015 Shangcheng Co., Anhui, China KY800607  KC147811 
90 Tylototriton dabienicus lineage 1 HNNU 1004-024 Shangcheng Co., Anhui, China KY800608 КС147812 
91 Tylototriton dabienicus lineage 1 HNNU 1004-026 Shangcheng Co., Anhui, China KY800609 КҮ800869 
92 it Misi У СІВ WH10001 Wufeng Со., Hubei, China KY800600 KY800863 
93 i id d EOS UE CIB WH10002 Wufeng Co., Hubei, China KY800601 — KY800864 
94 I pie аах CIB WH10003 Wufeng Co., Hubei, China KY800602 — KY800865 
95 ee ак СІВ WH10007 Wufeng Co., Hubei, China KY800603 KY800866 
96 " ` osi СІВ Wg20090730001 Libo Co., Guizhou, China KY800575 KY800842 
97 n жө TE СІВ Wg20090730002 Libo Co., Guizhou, China KY800576  KY800843 
98 t еа а ы СІВ Wg20090730003 Libo Co., Guizhou, China KY800577 KY800844 
99 " E dion E CIB Wg20090730005 Libo Co., Guizhou, China KY800578 КҮ800845 
100 " nid wees СІВ Wg20090730004 Libo Co., Guizhou, China KY800580 KY800846 
101 He plas wenxtanensis СІВ WG200600019 ^ Suiyang Co., Zunyi, Guizhou, China KY800544 KY800817 
102 een i CIB WG20060007 Suiyang Co., Zunyi, Guizhou, China KY800545 Кү800818 
103 in Pin iE CIB WG20090601002  Suiyang Co., Zunyi, Guizhou, China KY800573 — KY800840 
104 He — а СІВ WG20090601001  Suiyang Co., Zunyi, Guizhou, China KY800574 KY800841 
105 oe аа СІВ 20090527 Wenxian Co., Gansu, China KY800579 КС147813 
106 n “н МЕИ СІВ 2010123101 Pingwu Со., Gansu, China KY800604 КҮ800867 
107 io PEERS CIB 2010123102 Pingwu Co., Gansu, China KY800605 Кү800868 
108 Иол Welixtanensis CIB 20070639 Qingchuan Co., Sichuan, China KY800542 KY800815 


lineage 1 
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Continued 





ae Species name Voucher number Locality 165 rRNA ND2 

109 ү, 2i Pid ы CIB 20070638 Qingchuan Co., Sichuan, China KY800543 Кү800816 
110 |, ы we d СІВ 20080002 Yunyang Co., Chongqing, China KY800540 KY800813 
111 " iin, тееп CIB 20080003 Yunyang Co., Chongqing, China KY800541 KY800814 
112 ко i CIB 20081201 Yunyang Co., Chongqing, China KY800567 KY800835 
113 шек iid CIB 20081202 Yunyang Co., Chongqing, China KY800568 Кү800836 
114 лен wenxianensis CIB WA20090601 Wangcang Co., Sichuan, China KY800571 КҮ800838 
115 " н аи СІВ WA20090602 Wangcang Со., Sichuan, China KY800572 KY800839 

Outgroup: 

116 Echinotriton chinhaiensis CIB ZHJY1 Zhenhai Co., Zhejiang, China KY800627 KY800891 
117 Echinotriton chinhaiensis CIB ZHJY2 Zhenhai Co., Zhejiang, China KY800628 KY800892 
118 Echinotriton andersoni MVZ 232187 Tokunoshima, Kagoshima, Japan EU880314 EU880314 
119 Pleurodeles waltl MVZ 231894 Cadiz, Andalusia, Spain EU880330 EU880330 


For sampling localities see Figure 1. Institutional abbreviations: CAS: California Academy of Sciences, Department of Herpetology, USA; CIB: 
Chengdu Institute of Biology, Chinese Academy of Sciences, China; CSUFT: Central South University of Forestry and Technology, China; FMNH: 
Field Museum, USA; HNNU: Henan Normal University, China; IEBR: Institute of Ecology and Biological Resources, Vietnam; KIZ: Kunming Institute 
of Zoology, Chinese Academy of Sciences, China; KUHE: Graduate School of Human and Environmental Studies of Kyoto University, Japan; 
LSUHC: La Sierra University, Herpetological Collection, USA; MVZ: Museum of Vertebrate Zoology, University of California, USA; TBU: Thuan Chau 
University, Vietnam; VNMN: Vietnam National Museum of Nature, Vietnam; VNUH: Vietnam National University Hanoi, Vietnam; ZDUM: Zoological 
Department, University of Mandalay, Myanmar; ZISP: Zoological Institute, Russian Academy of Sciences, Russia; ZMMU: Zoological Museum of 


Lomonosov Moscow State University, Russia. Co.: County. —: Not available. 


Phylogenetic analyses 

Sequences of partial fragments of ND2 and 16S rRNA mtDNA 
for 119 Salamandridae specimens, including 115 representatives 
of Tylototriton (26 species) and four sequences of outgroup 
members of Salamandridae (Echinotriton and Pleurodeles) 
were included in the final alignment with a total length of up to 
1 665 bp. Information on voucher specimens and GenBank 
accession Nos. used in phylogenetic analyses is summarized 
in Table 1. Nucleotide sequences were initially aligned in 
MAFFT v.6 (Katoh et al., 2002) with default parameters, and 
then checked by eye in BioEdit 7.0.5.2 (Hall, 1999) and 
slightly adjusted. 

The dataset was divided into four partitions: three codon- 
partitions for the ND2 gene and a single partition for 16S 
rRNA, with the optimal evolutionary models for each estimated 
using MODELTEST v. 3.06 (Posada & Crandall, 1998). 
According to the Akaike information criterion (AIC), the TVM+ 
G model was the best fit for the 16S rRNA partition; for the 
ND2 gene, however, the HKY * G model was considered the 
best fit for the first and second codon partitions, whereas the 
J2+G model was selected as the best fit for the third codon 
partition. Mean uncorrected genetic distances (P-distances) 
between sequences were determined with MEGA 7.0 (Kumar 


et al., 2016). 

The matrilineal genealogy was inferred using Bayesian 
inference (BI) and maximum likelihood (ML) algorithms. BI 
analyses were conducted in MrBayes v3.1.2 (Huelsenbeck & 
Ronquist, 2001; Ronquist & Huelsenbeck, 2003). Metropolis- 
coupled Markov chain Monte Carlo (MCMCMC) analyses 
were run with one cold chain and three heated chains for 
twenty million generations and sampled every 2 000 
generations. Five independent MCMCMC runs were 
performed and 1 000 trees were discarded as burn-in. We 
checked the convergence of the runs and that the effective 
sample sizes (ESS) were all above 200 by exploring the 
likelihood plots using TRACER v1.6 (Rambaut et al., 2014). 
Confidence in tree topology was tested by posterior probability 
(PP) for the BI trees (Huelsenbeck & Ronquist, 2001). Nodes 
with PP values over 0.95 were a-priori regarded as sufficiently 
resolved, those between 0.95 and 0.90 were regarded as 
tendencies, and values below 0.90 were considered to be not 
supported. 

We conducted ML analyses using the RAxML web server 
(http://embnet. vital-it. ch / raxml-bb/; Stamatakis et al., 2008) 
and searched ML trees using the gamma model of rate 
heterogeneity option. Confidence in node topology was tested 
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by non-parametric bootstrapping with 1 000 replicates (ML 
BS, see Felsenstein, 1985). We a-priori regarded tree nodes 
with bootstrap (ML BS) values of 70% or greater and BI PP 
values over 0.95 as sufficiently resolved; ML BS values 
between 70% and 50% (BI PP between 0.95 and 0.90) were 
treated as tendencies, and nodes with ML BS values below 
50% (BI PP below 0.90) were regarded as unresolved 
(Felsenstein, 2004; Huelsenbeck & Hillis, 1993). 


RESULTS 


Phylogenetic analyses 

Sequences and statistics: The final alignment of the ND2 
gene contained 1 157 aligned characters, including 712 
conserved sites and 445 variable sites, of which 405 were 
parsimony-informative. The transition-transversion bias (R) was 
estimated to be 4.56 (all data for ingroup only). Nucleotide 
frequencies were 37.5% (A), 23.7% (T), 28.3% (C), and 
10.5% (G). The final alignment of the 16S rRNA gene 
contained 508 aligned characters, including 424 conserved 
sites and 82 variable sites, of which 69 were parsimony- 
informative. The transition-transversion bias (R) was estimated 
to be 5.84 (all data for ingroup only). Nucleotide frequencies 
were 36.896 (A), 24.996 (T), 20.396 (C), and 18.096 (G). 


Position of Tylototriton sp. in matrilineal genealogy: В! 
and ML phylogenetic analyses resulted in essentially similar 
topologies (Figure 2). In general, the topology of the mtDNA- 
based matrilineal genealogy was consistent with the phylogeny 
of Tylototriton presented by Wang et al. (2018), suggesting 
that the genus is divided into five clades (1—5) grouped into 
two major reciprocally monophyletic groups (Figure 2): 

(1) The first group joined two clades with sister relationships: 
clade 1 (including T. verrucosus (type species of Tylototriton s. 
str. Anderson, 1871), T. anguliceps, T. himalayanus, T. 
kweichowensis (type species of Qiantriton Fei, Ye & Jiang, 
2012), T. ngarsuensis, T. podichthys, T. pulcherrimus, T. 
shanjing, T. shanorum, T. uyenoi, T. yangi,and Tylototriton sp. 
1 from the western part of the Kachin and Sagaing States of 
Myanmar and the newly discovered population of Tylototriton 
sp. nov. from the Indawgyi Lake area in the southern part of 
Kachin State) and clade 2 (including T. pseudoverrucosus and 
T. taliangensis, the latter being the type species of 
Liangshantriton Fei, Ye & Jiang, 2012). 

(2) Thesecond group joined clade 3 (including T. broadoridgus, 
T. dabienicus, T. liuyangensis, T. lizhenchangi, and T. 
wenxianensis), clade 4 (including T. panhai and T. 
vietnamensis), and clade 5 (including T. asperrimus (type 
species of Yaotriton Dubois & Raffaélli, 2009), T. hainanensis, 
T. notialis, and T. ziegleri); the topological relationships 
between these three clades are essentially unresolved. 

In accordance with the results of Wang et al. (2018) our 
analysis indicated deep phylogenetic structuring and paraphyly 
of T. asperrimus (consisting of two non-monophyletic lineages), 
T. wenxianensis (consisting of three lineages, not forming a 
monophyly), and Т. dabienicus (consisting of two non- 
monophyletic lineages), suggesting that taxonomy of this 
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group is incomplete and further taxonomic and phylogenetic 
research is required. 

The newly discovered population of Tylototriton sp. nov. 
from the Indawgyi Lake area belongs to clade 1 (Figure 2), 
which occurs in western and northern Indochina, Himalaya, 
and Yunnan Province of China (Figure 1), and is grouped with 
T. himalayanus from Nepal, though not with significant node 
support (0.85/57, hereafter given for ВІ PP/ML BS, 
respectively). Tylototriton sp. nov. from the Indawgyi Lake 
area and T. himalayanus form a well-supported monophyletic 
group with Tylototriton species from the Shan Plateau in 
Myanmar (1.0/96), whereas populations of Tylototriton sp. 1 
from the Sagaing and Kachin states in northern Myanmar, 
reported by Grismer et al. (2018a), are clustered with other 
members of Tylototriton clade 1 from Yunnan and northern 
Indochina and belong to the 7. verrucosus species complex 
(Figure 2). Our analyses suggest sister relationships of closely 
related T. verrucosus and T. shanjing but only with support 
from BI (1.0/-), suggesting that the taxonomic status of T. 
shanjing might need to be reconsidered. Finally, in our tree, 
the recently described Myanmar species T. ngarsuensis is 
nested within differentiation of T. shanorum, rendering the 
latter paraphyletic. 


Sequence divergence: The uncorrected P-distances among 
and within the studied mtDNA fragments for the examined 
Tylototriton species are shown in Table 2 (data for ingroup 
only). The interspecific uncorrected genetic P-distances between 
the Tylototriton sp. from Kachin State of Myanmar and other 
congeners varied from 5.396 (between Tylototriton sp. and its 
sister species T. himalayanus) to 14.6% (between Tylototriton 
sp. and T. lizhenchangi) for the ND2 gene; and from 2.496 
(between Tylototriton sp. and its sister species T. himalayanus) 
to 5.996 (between Tylototriton sp. and T. liuyangensis and T. 
dabienicus, lineage 1) for the 16S rRNA gene (Table 2). This 
degree of pairwise divergence is quite high, notably greater 
than the genetic divergence observed between many 
recognized species of Tylototriton (see Table 2 and Grismer et 
al., 2018a; Wang et al., 2018). 


Taxonomy 

Our mtDNA genealogy analyses based on the ND2 and 16S 
rRNA genes indicated that the newly discovered population of 
Tylototriton sp. nov. from the Indawgyi Lake area belongs to 
clade 1 of the subgenus Tylototriton s. str. and is clustered 
with two other species of the genus known from Myanmar, T. 
ngarsuensis and T. shanorum, and with T. himalayanus from 
Nepal (Figure 2). The lineage of Tylototriton sp. nov. from 
Indawgyi Lake is clearly distinct and notably divergent from all 
other congeners with the uncorrected genetic distance for 
interspecific comparisons exceeding P=5.3% in the ND2 gene 
and P=2.4% in the 16S rRNA gene. The observed differences 
in mtDNA sequences are congruent with evidence from 
diagnostic morphological characters (see "Comparisons"). 
These results support our hypothesis that the newly discovered 
population of Tylototriton sp. nov. from Indawgyi Lake represents 
a previously unknown species, which we describe herein. 
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Figure 2 Bayesian inference consensus tree of genus Tylototriton derived from analysis of 1 157 bp ND2 and 508 bp 16S rRNA gene 


fragments 
Specimen numbers (1—115), corresponding voucher specimen information, and GenBank accession Nos. are given in Table 1. Color denotes five 


major clades revealed within genus Tylototriton (clades 1—5). Information on type species of subgeneric-level taxa Tylototriton s. str., Qiantriton, 
Liangshantriton, and Yaotriton is provided. Number at tree nodes corresponds to BI PP/ML BS support values, respectively. Photo by Nikolay A. 


Poyarkov. 
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Tylototriton kachinorum sp. nov. 
Tables 3-4; Figures 3-6. 


Holotype: ZMMU A5953 (field number NAP-08318), adult 
male from a swamp in a forest clearing surrounded by montane 
evergreen tropical forest, Ingyin Taung Mountain, Indawgyi 
Lake area, Mohnyin Township, Kachin State, Myanmar 
(approximate coordinates N25.09°, E96.28°; elevation 1 000 
m а.5.1.), collected on July 18, 2018 at 2100 h by Than Zaw, 
Paw Lay, Parinya Pawangkhanant, Vladislav A. Gorin, and 
Nikolay A. Poyarkov. 


Paratypes: ZMMU A5954 (field number NAP-08320), ZISP 
13721 (field number NAP-08324), ZDUM-0101-0105 (field 
numbers NAP-08325, NAP-08322, NAP-08319, NAP-08326 
and NAP-08317, respectively), seven adult males from the 
same locality and with the same collection information as the 
holotype; and ZMMU А5955 – А5956 (field numbers NAP- 
08323 and NAP-08321, respectively), two adult females from 
the same locality and with the same collection information as 
the holotype. 


Referred specimens: ZMMU A5957 (field number NAP-08305), 
a larva (Grosse (2013) stage 40) from the same locality and 
with the same collection information as the holotype. 


Diagnosis: The new species is assigned to the genus 
Tylototriton based on molecular data and by the following 
combination of morphological attributes: (1) presence of dorsal 
granules, (2) dorsolateral bony ridges on head, (3) presence 
of dorsolateral series of knob-like warts (rib nodules); and 
(4) absence of quadrate spine ( Figure 2). Tylototriton 
kachinorum sp. nov. is distinguished from all other congeners 
by a combination of the following morphological attributes: 
(1) medium body size, adult SVL 62.3-74.1 mm in males, 72.5— 
84.8 mm in females; (2) tail thin and long, longer than body in 
both sexes, lacking lateral grooves; (3) skin rough with fine 
granules; (4) snout truncate in dorsal view; (5) supratemporal 
bony ridges on head wide, protruding, beginning at anterior 
corner of orbit; (6) sagittal ridge on head very weak, almost 
indistinct; (7) limbs long and thin, tips of forelimb and hindlimb 
broadly overlapping when adpressed along body; (8) vertebral 
ridge distinct, wide, non-segmented; (9) rib nodules weakly 
distinct, 13 — 14 along each side of body; (10) background 
coloration brown to dark-brown; (11) labial regions, parotoids, 
rib nodules, whole limbs, vent, ventral tail ridge with dull 
orange-brown to yellowish-brown markings. 

The new species is also markedly distinct from all congeners 
for which comparable sequences are available of ND2 (P2 
5.3%) and 16S rRNA (P22.496) mitochondrial DNA genes. 


Description of holotype: A medium-sized specimen in a 
good state of preservation (Figures 3—4). 


Head: Head longer than wide (HW/HL ratio 89.396) (Figure 3C), 
head wider than body; angularly hexagonal in shape in dorsal 
view, slightly depressed, gently sloping in profile (Figure 3E); 
snout comparatively long, three times longer than eye (UEW/ 


SL ratio 33.996), sharply truncate in dorsal view (Figure 3C), 
slightly rounded in lateral view (Figure 3E), slightly projecting 
beyond lower jaw; nostrils on anterior margin of snout located 
notably closer to snout tip than to eye (NSD/ON ratio 58.8%), 
facing anterolaterally, not visible from dorsal view; labial folds 
absent; tongue oval, attached to anterior floor of mouth but 
free posteriorly and laterally; vomerine teeth arranged in 
inverted V-shaped almost straight series, converging and 
narrow anteriorly, gradually widening posteriorly, notably 
longer than wide (VTW/VTL ratio 77.796), anteriorly reaching 
beyond level of choanae but not in contact with them, vomerine 
teeth 95 (47/48 in right and left branches, respectively), upper 
jaw teeth 93, and lower jaw teeth 110; parotoids distinct, 
comparatively large, crescent-shaped, notably projecting 
posteriorly (Figure 3E); dorsolateral supratemporal bony 
ridges on head wide, notably protruding, from anterior corner 
of orbit to anterior end of parotoid, forming medially recurved 
projection on its posterior end (Figure 3C); sagittal bony ridge 
on head very weak, almost indistinct (Figure 3C); gular fold 
present (Figure 3D). 


Body: Body habitus comparatively slender (Figure 3A); costal 
folds absent; vertebral middorsal ridge wide, non-segmented, 
running from occiput region to anterior one fifth of tail length, 
separated from sagittal head ridge on head with wide gap 
(Figure 3C); rib nodules weakly distinct, small, forming knob- 
like glandular warts, arranged in two longitudinal lines on 
dorsolateral surfaces of dorsum, 14 on both sides of body 
from area posterior to axilla to level of posterior vent margin 
(base of tail) (Figure 3A); rib nodules almost of same size, 
rounded, those in posterior third of dorsum slightly oval- 
shaped, those on sacral area notably elongated, decreasing in 
size posteriorly on sacrum and tail basis. Limbs: Limbs 
comparatively long, slender (Figure 3A); forelimbs slightly 
shorter than hindlimbs; relative length of forelimb FLL/SVL 
ratio 26.2%, relative length of hindlimb ratio 28.0%; fore- and 
hindlimbs largely overlapping when adpressed towards each 
other along sides of body; fingers and toes well developed 
(Figure 3F-1), free of webbing; fingers four, comparative finger 
lengths: 1FL<4FL<2FL<3FL; toes five, comparative toe 
lengths: 1TL<5TL<2TL<3TL<4TL. Tail: Tail very long, notably 
exceeding body length (TAL/SVL ratio 114.9%); tail laterally 
compressed along entire length, tapering posteriorly, lateral 
grooves on tail absent; dorsal tail fin starting at anterior one 
fifth of tail length, more distinct posteriorly, with maximal tail 
height at posterior two thirds of tail length, dorsal tail fin 
slightly serrated; ventral tail fin smooth; tail tip pointed. Skin 
texture and skin glands: Skin rough, small granules present 
on dorsal surfaces of head and dorsum (Figure 3A,C), lateral 
sides of body and tail; on ventral surface granules become 
smaller, arranged in transverse striations (Figure 3B); small, 
sparse granules regularly arranged on throat (Figure 3D); 
head ridges with rough surface; skin on volar and plantar 
surfaces of hands and feet with tiny grooves forming 
reticulated pattern; metacarpal or metatarsal tubercles absent. 
Cloacal region notably swollen, vent as longitudinal slit (Figure 
3J), vent edges with numerous transverse folds. 
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Figure 3 Holotype of Tylototriton kachinorum sp. nov. (ZMMU A5953, male) in life 
A: Dorsal view; B: Ventral view; C: Head, dorsal view; D: Head, ventral view; E: Head, lateral view; F: Opisthenar view of right hand; G: Volar view of 
right hand; H: Opisthenar view of right foot; l: Plantar view of right foot; J: Ventral view of cloacal area. Photos by Nikolay A. Poyarkov. 
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Color of holotype in life: Dorsal ground color of dorsal surfaces 
of head and trunk dark brown(Figure4); dorsal surfaces of limbs 
and lateral surfaces of tail light yellowish-brown (Figure 3A); 
iris brown with tiny bronze speckles along outer margins 
(Figure 3E); gular region, belly, and ventral surfaces of limbs 


light yellowish-gray (Figure 3B); anterior parts of head and 
parotoids light orange-brown; rib nodules and vertebral ridge 
yellowish-brown to orange-brown, barely discernable from dark 
brown trunk coloration; upper and lower lips, posteriormost 
corners of parotoids, palms, and soles light-orange to yellowish. 





Figure 4 Holotype of Tylototriton kachinorum sp. nov. (ZMMU A5953, male) in situ (Photo by Nikolay A. Poyarkov and Than Zaw) 


Color of holotype in preservative: After preservation in ethanol 
for six months, coloration pattern of holotype resembles that 
observed in life, however yellowish and orange tints faded to 
light brownish-gray. 


Measurements and counts of holotype: Morphometric 
characters (all in mm): SVL 67.2; RHL 18.1; RHW 16.2; 
RMXHW 16.6; RIND 5.4; RAGD 34.1; RTRL 50.4; TAL 77.3; 
RVL 6.3; RFLL 26.2; RHLL 28.0; RVTW 6.0; RVTL 7.7; RLJL 
13.5; RSL 6.8; RIOD 7.2; RUEW 2.3; RUEL 4.4; ROL 3.3; 
RBTAW 9.0; RMTAW 3.2; RMXTAH 10.8; RMTAH 7.0; RON 
4.6; ICD 9.6; CW 11.5; NSD 2.7; 1FL 2.9; 2FL 4.8; 3FL 5.5; 
4FL 3.0; 1TL 2.5; 2TL 5.4; 3TL 7.0; 4TL 7.7; 5TL 3.6. Meristic 
characters: UJTN 93; LUTN 110; VTN 47/48 (right/left). 


Variation: All individuals in the type series were generally similar 
in morphology and agreed with the holotype description in 
body proportions and coloration; variation of morphometric 
characters within the type series is shown in Table 3. The 
variation of dorsal coloration in seven male and two female 
paratypes in life is presented in Figure 5. In general, males 
had more robust and slender bodies than females. The males 
were notably smaller in body size (SVL 62.3-74.1 mm, mean 
68.6 € 2.9 mm) than the two females (SVL 72.5 – 84.8 mm) 
(Table 3). Body of the largest female (ZMMU A5956) was 
notably swollen and was wider than head width (Figure 5). 


Male paratypes ZDUM-0105 and ZISP 13721 had notably shorter 
tails than other type specimens (Table 4) due to regeneration 
of tail tip after damage (Figure 5). Coloration within the type 
series slightly varied, from specimens lighter than the holotype, 
which appeared light-brown to yellowish-orange (ZDUM-0102), 
to darker specimens, which appeared dark-brown with duller 
orange-brown light markings (male ZISP 13721 and female 
ZMMU А5956) (Figure 5). 


Eggs and clutch: The clutch size is unknown. The diameter 
of ripe eggs in the ovaries of a female paratype (ZMMU A5956) 
ranged from 1.6 to 1.7 mm (n=5, mean=1.7 mm). The animal 
pole was dark brown and the remaining area of ova was cream. 


Larval morphology: Description of larval morphology is 
based on a single larval specimen (ZMMU A5957, Grosse 
(2013) stage 40) (see Referred specimens for details). Details 
of larval morphology are presented in Figure 6. 


Larval measurements (all in mm): SVL 10.5; HL 3.8; HW 
3.4; OL 0.8; AGD 5.9; TAL 9.7; FLL 3.6; HLL 2.8; MXTAH 2.6. 


Larval external morphology: Body elongated, higher than 
wide. Head large, trapezoidal in shape, wide and slightly 
depressed with short, flattened snout, comprising 3796 of 
snout-vent length (SVL), gently sloping in lateral view, two 
times wider than body in dorsal view. Snout truncate in dorsal 
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ZDUM-0104 ZDUM-0103 ZISP 13721 ZDUM-0102 ZDUM-0105 


ZDUM-0101 ZMMU A5954 ZMMU A5955 ZMMU A5956 


Figure 5 Variation of dorsal coloration in paratypes of Tylototriton kachinorum sp. nov. 
Scale bar: 10 mm. Photos by Nikolay A. Poyarkov. 
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Table 3 Measurements of type series of Tylototriton kachinorum sp. nov. (all in mm) 








Specimen ZMMU ZMMU ZDUM- ZDUM- ZDUM- ZDUM- ZDUM- ZISP ZMMU ZMMU 

А5953 А5954 0101 0102 0103 0104 0105* 13721* A5955 — A5956 — F 
ow Holotype Paratype Paratype Paratype Paratype Paratype Paratype Paratype Paratype Paratype 
Sex M M M M M M M M way (Min-Max) F F jos 
SVL 67.2 74.1 74.0 68.3 66.9 69.3 62.3 66.3 68.6+2.9 (62.3-74.1) 72.5 84.8 78.7+6.2 
HL 18.1 20.7 20.7 18.6 19.6 19.7 16.9 17.3 18.9+1.2 (16.9-20.7) 18.1 20.3 19.241.1 
HW 16.2 19.2 18.5 17.0 16.7 VT 15.8 16.4 17.2+1.0 (15.8-19.2) 18.4 21.4 19.9+1.5 
MXHW 16.6 19.6 18.6 17:5 17.6 18.4 16.6 16.7 17.7+0.9 (16.6—19.6) 19.4 23.0 21.2+1.8 
IND 5.4 6.4 5.7. 5.4 5.2 5.6 4.2 4.8 5.340.4 (4.2—6.4) 5.8 6.8 6.3.+0.5 
AGD 34.1 36.6 36.6 34.0 33.1 33.9 30.8 33.2 34.0+1.3 (30.8-36.6) 40.5 48.2 44.4+3.9 
TRL 50.4 56.2 55.1 51.3 50.4 52.9 47.1 49.4 51.6+2.3 (47.1-56.2) 58.4 69.0 63.7+5.3 
TAL 77.3 89.6 96.3 79.2 79.2 85.2 57.0 67.8 78.948.7 (57.0-96.3) 82.7 85.2 84.0+1.3 
VL 6.3 7.4 9.7 6.8 6.2 7.4 5.7 6.5 7.040.9 (5.7—9.7) 4.9 4.8 4.90.1 
FLL 26.2 29.3 29.1 26.5 26.1 28.2 24.1 26.4 27.0+1.4 (24.1—29.3) 26.1 29.9 28.0+1.9 
HLL 28.0 31.0 31.6 28.0 27.3 29.7 26.6 27.5 28.7+1.5 (26.6-31.6) 29.0 32.2 30.6+1.6 
VTW 6.0 6.7 6.6 6.2 6.1 6.3 5.6 5.9 6.2*0.3 (5.6—6.7) 6.4 7.4 6.9+0.5 
VTL TA 9.1 9.4 8.8 8.9 9.2 8.4 8.8 8.840.4 (7.7—9.4) 8.4 9.7 9.1+0.6 
LJL 13.5 16.0 14.1 13.9 13.7 14.4 12.3 13.4 13.9+0.7 (12.3-16.0) 14.5 16.9 15.741.2 
SL 6.8 7.6 8.1 7.3 7.0 7.5 6.2 6.5 7.140.5 (6.2-8.1) 7.3 8.7 8.0+0.7 
IOD 7.2 8.4 8.1 7.3 8.0 7.9 7.2 7.3 7.740.4 (7.2—8.4) CH 8.7 8.2+0.5 
UEW 2:3 2.5 2.6 2.4 2.4 2.4 2.2 2.2 2.4+0.1 (2.2-2.6) 2.3 2.8 2.6+0.2 
UEL 4.4 5.0 4.9 4.6 4.4 4.6 4.1 4.4 4.540.2 (4.1-5.0) 4.3 5.3 4.8+0.5 
OL 3.3 3.6 3.6 3.3 3.2 3.4 3.1 3.2 3.340.1 (3.1-3.6) 3.6 4.2 3.9+0.3 
BTAW 9.0 9.7 9.5 10.0 9.7 9.6 8.8 9.6 9.5*0.3 (8.8-10.0) 8.5 9.8 9.2+0.7 
MTAW 3.2 3.3 3.9 3.4 3.6 3.7 3.3 3.6 3.5+0.2 (3.2—3.9) 3.8 4.2 4.0+0.2 
МХТАН 10.8 11.8 10.8 10.7 10.1 10.2 9.0 10.9 10.5+0.6 (9.0-11.8) 11.6 11.9 11.7+0.2 
МТАН 7.0 9.0 7.3 7.5 7.8 CL 7.7 8.4 7.840.5 (7.0—9.0) 9.3 9.7 9.5+0.2 
ON 4.6 5.1 5:1 4.7 4.7 4.7 4.2 4.7 4.740.2 (42—51) 5.0 5.2 5.1+0.1 


Character abbreviations: SVL: Snout-vent length; HL: Head length; HW: Head width; MXHW: Maximum head width; IND: Internarial distance; AGD: 
Axilla-groin distance; TRL: Trunk length; TAL: Tail length; VL: Vent length; FLL: Forelimb length; HLL: Hindlimb length; VTW: Vomerine tooth series 
width; VTL: Vomerine tooth series length; LUL: Lower jaw length; SL: Snout length; IOD: Interorbital distance; UEW: Upper eyelid width; UEL: Upper 
eyelid length; OL: Orbit length; BTAW: Basal tail width; MTAW: Tail width at mid-level of tail; MXTAH: Maximum tail height; MTAH: Tail height at mid- 
level of tail; ON: Orbitonarial distance. M: Male. F: Female. For other abbreviations see “Materials and Methods”. Asterisk (*) denotes damaged and 


regenerated tail. 


view (Figure 6B), rounded in lateral view (Figure 6A). Tail 
subequal to body length comprising 93% of SVL; myotomes 
on body and tail not discernable in lateral view. Nostrils 
rounded, small, oriented anterolaterally, located much closer 
to snout tip than to eye. Eyes large, rounded, with lateral 
orientation but still visible in dorsal view (Figure 6B). Limbs 
thin, forelimbs longer than hindlimbs, HLL/FLL ratio 79.1%. 
Forelimbs with four well-developed elongated fingers; forelimb 
turned, palm facing ventrally; relative finger lengths: 4FL<3FL< 
1FL«2FL. Hindlimbs with knee joint already formed and four 


well-developed toes, fifth toe as nub; relative toe lengths: 5TL< 
ATL<3TL<1TL<2TL. Orbit diameter (OL) 7.7% of SVL. Short 
longitudinal slit for vent. Height of tail musculature at highest 
portion comprises 50%—60% of tail height. Maximum height of 
dorsal tail fin 4596—5596 of maximum tail height. Ventral tail fin 
two times lower than dorsal tail fin. Ventral tail fin roughly at 
level of vent, dorsal fin lower than head, at level of axilla and 
reaching maximum height mid tail. Tail tip sharply pointed 
(Figure 6A). Skin completely smooth; lateral line organs visible 
on ventral side of head; mouth open with well-developed 
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Figure 6 Lateral (A), dorsal (B), and ventral (C) views of larval specimen (ZMMU A5957; Grosse (2013) stage 40) of Tylototriton 


kachinorum sp. nov. in life 
Scale bar: 3 mm. Photos by Nikolay A. Poyarkov. 
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teeth; no remains of yolk (Figure 6C); gills well-developed, 
much higher than body, with fimbriae clearly visible. 


Larval coloration in life: In life body background color ochre 
to golden (Figure 6A, B), ventral surfaces pinkish, translucent 
(Figure 6C). Body, tail, and head pigmented dorsally: tail 
almost uniform purple-gray with rare golden speckles, 
pigmentation forms tortoise-shell golden—dark-gray pattern on 
body and dorsal fin, head pigmentation less dense, reduced to 
grayish spots and dots. Few dark spots on ventral fin and 
limbs. Eyes, except for pupil, fully pigmented, iris golden 
(Figure 6A). 


Position in mtDNA genealogy and sequence divergence: 
According to our mtDNA data, Tylototriton kachinorum sp. 
nov. belongs to clade 1 of the subgenus Tylototriton s. str. 
(Figure 2) and is grouped with Tylototriton species from the 
Shan Plateau of Myanmar (T. shanorum, T. ngarsuensis) and 
Himalaya (T. himalayanus). Uncorrected genetic P-distances 
between Tylototriton kachinorum sp. nov. 16S rRNA 
sequences, and all homologous sequences of congeners 
available included in our analyses varied from 5.3% (with 
sister species T. himalayanus) to 14.6% (with T. lizhenchangi) 
(Table 2). 


Distribution and biogeography: To date, Tylototriton 
kachinorum sp. nov. is known only from a single locality on 
the slopes of Ingyin Taung Mountain, Kachin State, Myanmar 
(Figure 1) at elevations from 900 to 1 050 m a.s.l. The Ingyin 
Taung Mountain belongs to the southernmost part of the 
Kachin Hills — a heavily forested group of highlands in the 
extreme northeastern area of Kachin State, consisting of a 
series of mountain ranges running mostly in a north-to-south 


direction. The actual distribution of Tylototriton kachinorum sp. 
nov. may be wider: it is anticipated that the new species 
occurs in the montane forests of adjacent mountains 
surrounding the largest inland lake of Kachin State, Indawgyi 
Lake, and possibly further northwards along the Kumon Bum 
subrange of the Kachin Hills. 


Natural history notes: Our knowledge on the biology of 
Tylototriton kachinorum sp. nov. is scarce. Adult animals were 
encountered at night after 1900 h in flooded areas of shallow 
slow-moving streams and artificial ponds in forest clearings, 
which were used by local Kachin farmers as a watering place 
for cattle (Figure 7). Surrounding areas were covered by 
secondary bamboo forest and primary mixed evergreen 
tropical forest. Adult male newts were observed slowly moving 
along the clay bottom in clear water 20—40 cm deep; both 
females and the larval specimen were collected in deeper 
areas (60—100 cm deep) from dense water vegetation using a 
dip-net. Courtship behavior of male newts was observed in 
July. Local Kachin farmers reported that they often find adult 
newts walking in dense vegetation far from waterbodies, 
especially after rain; they are also often encountered in wells 
they construct. The new species is known to local Kachin 
people as “Lan Yan" literally meaning "water agama” in their 
native language. 

Other species of amphibians recorded syntopically with the 
new species at the type locality include Microhyla heymonsi 
Vogt, Microhyla mukhlesuri Hasan, Islam, Kuramoto, 
Kurabayashi & Sumida, Microhyla butleri Boulenger, 


Limnonectes limborgi (Sclater), Limnonectes sp., Fejervarya 
sp., Kurixalus sp., Feihyla vittata (Boulenger), Rhacophorus 
bipunctatus Ahl, Polypedates mutus (Smith), and Raorchestes 
parvulus (Boulenger). 





Figure 7 Breeding habitats of Tylototriton kachinorum sp. nov. at type locality 
Temporary swamps in forest clearings within montane tropical forest on Ingyin Taung Mountain, Indawgyi Lake area, Kachin State, Myanmar, where 
larva and adult specimens were collected. Photos by Parinya Pawangkhanant. 
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Comparisons: According to phylogenetic analyses, Tylototriton 
kachinorum sp. nov. falls into clade 1 of the subgenus 
Tylototriton s. str. and morphological comparisons with 
members of this group appear to be the most pertinent. The 
new species can be easily distinguished from members of the 
subgenus Yaotriton (clades 3-5 in Figure 2) by having light 
color markings on parotoids, lips, vertebral ridge, rib nodules, 
limbs, and ventral tail ridge (vs. dark body coloration except 
for palms and soles, vent region, and ventral ridge of tail in 
most members of the subgenus Yaotriton, with the exception 
of T. panhai). The new species can be further distinguished 
from T. panhai by having light color markings on entire limbs 
(vs. distinct light markings only on palms, soles, and fingers). 
Morphometric comparisons and morphological differences 
in several diagnostic characters between Tylototriton kachinorum 


sp. nov. and the closely related species of the subgenus 
Tylototriton are summarized in Table 4. In particular, the new 
species can be distinguished from T. taliangensis (in clade 2 
of the subgenus Tylototriton, sometimes regarded as a 
separate subgenus or full genus Liangshantriton; see Gong et 
al., 2018) by having distinct rib nodules, light markings on rib 
nodules, lips, and parotids (vs. lack of distinct rib nodules, 
mostly dark charcoal-black body coloration with light orange to 
red markings only on posterior part of parotoids, digits, palms, 
Soles, vent, and ventral tail ridge). From 7. pseudoverrucosus 
(clade 2) and T. kweichowensis (clade 1 of the subgenus 
Tylototriton, sometimes regarded as a separate subgenus 
Qiantriton), Tylototriton kachinorum sp. nov. can be 
distinguished by having isolated light markings on rib nodules 
(vs. connected markings forming light dorsolateral lines). 


Table 4 Morphological comparison between Tylototriton species found in Myanmar and adjacent territories 








Species Tylototriton kachinorum sp. nov. T. shanorum T. ngarsuensis 

M Е M F M F 
Character 8 2 1 2 2 1 
SVL 68.6+2.9 72.5-84.8 76.0 76.5-87.9 74.9-76.4 102.3 
RHL 27.6 23.9-24.9 22.4 24.3-25.7 24.0-26.0 21.5 
RHW 25.1 25.3-25.4 25.4 24.8-26.3 24.5-28.2 26.6 
RIND T 8.0 6.7 6.0-6.3 8.0-8.8 77 
RAGD 49.6 55.9-56.8 49.6 49.4-51.9 48.2-49.9 51.2 
RTRL 75.3 80.5-81.4 77.6 74.3-75.7 74.2-75.5 77.5 
RTAL 120.5 100.5-114.0 111.2 97.0-97.8 98.0-103.5 104.6 
RVL 10.2 5.6-6.8 9.3 3.4-3.5 10.7-12.3 8.0 
RFLL 39.3 35.3-36.0 34.7 31.7-32.5 39.7-39.9 35.5 
RHLL 41.9 37.9-39.9 37.2 35.6-37.0 41.9-47.1 38.7 
RMXHW 25.8 26.8-27.1 25.9 25.7-27.1 _ _ 
RVTW 9.0 8.7-8.8 6.1 6.1-7.5 = = 
RVTL 12.8 11.4-11.6 10.7 11.4-12.0 = = 
Snout Truncate Blunt to truncate Rounded to blunt 


Wide, protruding, begins at anterior cor- 


Supratemporal ridge ner of orbit 


Sagittal ridge Very weak, indistinct 


Surface of head ridges Rough 
Adpressed limbs Overlapping 
Vertebral ridge Wide, non-segmented 
Rib nodules Weakly distinct, 13-14 
Ground color Brown to dark brown 
Color of light markings Orange-brown to yellowish-brown 


Parotoids, rib nodules, palms, soles, 


Location of light markings vent, ventral tail ridge 


Lateral grooves on tail Absent 
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Wide, protruding, begins in loreal Wide, not protruding, begins pos- 


region terior to orbit 

Very weak Very weak, in males 
Rough Very rough 
Overlapping Overlapping 


Narrow, weakly segmented Wide, weakly segmented 


Distinct, 14 Distinct, 15 
Dark brown to black Nearly black 
Yellow to reddish brown Dark yellow 


Head, vertebral ridge, rib nodules, 
vent, whole limbs and tail 


Labial regions, palms and soles, 
vent, ventral tail ridge 


Absent Absent 


Continued 








Species T. himalayanus T. verrucosus T. shanjing 

M Е M M F M 
Character 32 13 3 32 13 3 
SVL 71.92+6.1 76.06+7.5 70.744.7 71.92+6.1 76.06+7.5 70.7+4.7 
RHL 24.5 25.9 24.3 24.5 25.9 24.3 
RHW 23.0 23.6 23.7 23.0 23.6 23.7 
RIND 8.4 8.2 7.0 8.4 8.2 7.0 
RAGD 53.6 51.6 50.7 53.6 51.6 50.7 
RTRL 79.9 77.1 75.7 79.9 77.1 (3:0 
RTAL 98.0 100.6 104.9 98.0 100.6 104.9 
RVL 11.8 12.5 7.0 11.8 12.5 7.0 
RFLL 36.2 36.1 34.3 36.2 36.1 34.3 
RHLL 37.5 38.1 37.3 37.5 38.1 37.3 
RMXHW Е _ 25.3 - - 25.3 
RVTW = - 7.7 = - 7.7 
RVTL = = 9.8 = = 9.8 
Snout Blunt Truncate Rounded 


Supratemporal ridge 
Sagittal ridge 

Surface of head ridges 
Adpressed limbs 
Vertebral ridge 

Rib nodules 

Ground color 


Color of light markings 


Location of light markings 


Very wide, protruding 
Weak, glandular 
Very rough 
Overlapping 
Non-segmented 
Large, prominent, 16 
Brown to dark brown 


Light brown 


Ventral surface of limbs, vent, ventral 


Narrow, steep 

Weak 

Smooth 

Overlapping 

Well segmented 
Weakly distinct, 14-15 
Blackish 

Orange 


Palms, soles, vent and ventral 


Narrow, steep 

Short, weak 

Smooth 

Overlapping 

Well segmented 
Distinct, 14 

Blackish 

Bright to dark orange 


Head, vertebral ridge, rib nodules, 








tail ridge ridge of tail vent, whole limbs and tail 
Lateral grooves on tail Very distinct Absent Weak 
Species T. uyenoi T. anguliceps T. podichthys 

M F M F M F 
Character 9 2 2 5 2 2 
SVL 68.1+3.8 69.3-78.3 61.1-62.5 70.6+3.4 56.5-60.2 73.4-78.3 
RHL 24.7 25.8-26.9 26.2-29.5 23.4 32.6-34.3 28.1—29.3 
RHW 25.0 23.1-24.0 22.7-23.4 22.7 26.4—28.0 24.8-25.9 
RIND 7.0 7.0-7.1 6.6-7.4 6:5 8.3-8.7 7.5-7.9 
RAGD 48.8 47.9-50.7 47.2-47.8 52.8 48.0—50.7 52.1-53.5 
RTRL 75.3 73.1-74.2 70.5-73.8 76.6 65.7-67.4 70.7-71.9 
RTAL 115.0 88.0-97.0 97.1-102.3 91.2 80.2-104.8 79.2-81.4 
RVL 7.4 1.7-1.9 T-2572 2.5 14.2-14.8 6.3-7.7 
RFLL 34.3 30.6-32.7 32.7-34.9 30.2 39.4-40.5 34.0-36.2 
RHLL 37.5 37.4-37.4 36.8-37.3 33.3 38.2-40.2 35.8-36.1 
RMXHW 25.9 23.6-25.0 23.7—24.1 23.3 = = 
RVTW 7.6 6.6-7.0 6.4-7.8 6.9 - _ 
RVTL 10.3 8.8-10.2 9.9-10.8 9.2 _ _ 
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Species T. uyenoi 


Continued 


T. anguliceps T. podichthys 





Snout Rounded 


Supratemporal ridge Narrow, steep 


Sagittal ridge Short, weak 
Surface of head ridges Smooth 
Adpressed limbs Overlapping 


Vertebral ridge Weakly segmented 


Rib nodules Prominent, 14—15 


Ground color Dark brown to black 


Color of light markings Dark-orange to red 


Head, vertebral ridge, rib nodules, vent, 


Location of light markings Whole limbs-and tail 


Lateral grooves on tail Absent 


Truncate Rounded 


Narrow, very steep Very wide, protruding 


Long, protruding Weak, glandular 


Rough Very rough 


Overlapping Touching 


Weakly segmented Non-segmented 


Prominent, 15 Large, prominent, 15—16 


Blackish Blackish 


Bright-orange Orange to dark-red 


Head, vertebral ridge, rib nodules, 
vent, dorsal surface limbs, ventral 
tail ridge 


Head, vertebral ridge, rib nodules, 
vent, whole limbs and tail 


Absent Absent 


Character abbreviations: SVL: Snout-vent length; RHL: Relative head length; RHW: Relative head width; RIND: Relative internarial distance; RAGD: 
Relative axilla-groin distance; RTRL: Relative trunk length; RTAL: Relative tail length; RVL: Relative vent length; RFLL: Relative forelimb length; 
RHLL: Relative hindlimb length; RMXHW: Relative maximum head width; RVTW: Relative vomerine tooth series width; RVTL: Relative vomerine 
tooth series length. M: Male. F: Female. For abbreviations see “Materials and Methods”. 


Tylototriton kachinorum sp. nov. can be distinguished from 
T. uyenoi, T. pulcherrimus, T. shanjing, and T. yangi by having 
dull orange-brown to yellowish-brown light markings (vs. much 
brighter orange to bright-yellow light markings). In particular, 
T. pulcherrimus has a series of bright-orange glandular spots 
located ventrolaterally and on flanks (vs. dark-brown coloration 
of flanks lacking light spots in new species); T. yangi has 
contrasting charcoal-black coloration of head and lips with 
only posteriormost part of parotoid colored bright orange, and 
no light ventral markings on body and tail (vs. all head dull 
orange-brown with slightly lighter lips and parotoids, light 
markings present on ventral tail ridge and vent in new 
species). Tylototriton kachinorum sp. nov. can be distinguished 
from T. uyenoi by relatively longer head in males (RHL 27.6 
vs. 24.7), greater internarial distance (RIND 7.7—8.0 vs. 7.0— 
7.1 for both sexes), and longer tail in both sexes (RTAL 120.5 
vs. 115.0 for males; 100.5— 114.0 vs. 88.0—97.0 for females) 
(Table 4). Males of the new species can be further diagnosed 
from males of T. shanjing by having comparatively wider head 
(RHW 25.1 vs. 22.2), longer tail (RTAL 120.5 vs. 104.4), 
greater internarial distance (RIND 7.7 vs. 7.1), comparatively 
wider (RVTW 9.0 vs. 6.7) and longer vomerine tooth series 
(RVTL 12.8 vs. 8.8) (Table 4). The new species can be further 
differentiated from T. shanjing by non-segmented vertebral 
ridge (vs. well-segmented) and brown to dark-brown background 
coloration of body (vs. blackish background coloration). 

Tylototriton kachinorum sp. nov. can be distinguished from 
T. verrucosus by having light ventral markings on body and tail 
(vs. no ventral markings on body and tail), comparatively 
longer head (RHL 27.6 vs. 24.3 in males, RHL 23.9—-24.9 vs 
21.7 in females), comparatively wider head in both sexes 
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(RHW 25.1—25.4 vs. 20.5-23.7), greater internarial distance in 
both sexes (RIND 7.7-8.0 vs. 6.2-7.0), longer tail in both 
sexes (RTAL 120.5 vs. 104.9 for males; 100.5-114.0 vs. 102.5 
for females) (Table 4), non-segmented vertebral ridge (vs. well- 
segmented), and brown to dark-brown background coloration 
of body (vs. blackish background coloration). 

The new species can be distinguished from T. podichthys 
by having comparatively shorter head in both sexes (RHL 23.9- 
27.6 vs 28.1—34.3), much longer tail in both sexes (RTAL 
100.5—120.5 vs. 79.2—104.8), truncate snout (vs. rounded), 
comparatively smoother skin on parotoids and dorsal surface 
of head, comparatively longer limbs (limbs widely overlap 
when adpressed to body in the new species vs. digit tips 
touch when limbs are adpressed to body in T. podichthys), 13— 
14 rib nodules (vs. 15—16 rib nodules), and duller coloration 
with orange-brown light markings and brown to dark-brown 
background (vs. orange to dark-red light markings and 
blackish background) (Table 4). 

Tylototriton kachinorum sp. nov. can be distinguished from 
T. anguliceps by having larger body size in both sexes (SVL 
68.6+2.9 mm in males and 72.5-84.8 mm in females of new 
species vs. 61.1—62.5 mm and 70.6+3.4 mm in T. anguliceps), 
comparatively wider head in both sexes (RHW 25.1—25.4 vs. 
22.7 — 23.4), greater internarial distance in both sexes (RIND 
7.7 —8.0 vs. 6.5-7.4), longer tail in both sexes (RTAL 100.5- 
120.5 vs. 91.2—- 102.3), and comparatively wider (RVTW 8.8— 
9.0 vs. 6.4-7.8) and longer vomerine tooth series (RVTL 11.4— 
12.8 vs. 9.2—10.8) (Table 4). The new species can be further 
diagnosed from T. anguliceps by having wide protruding 
supratemporal ridges (vs. narrow and steep supratemporal 
ridges), very small and almost indiscernible sagittal ridge (vs. 


long and notably protruding sagittal ridge), wide and non- 
segmented vertebral ridge (vs. narrow and weakly segmented 
vertebral ridge), less distinct 13-14 rib nodules (vs. rib 
nodules more distinct and protruding, not less than 15), and 
duller coloration with orange-brown light markings and brown 
to dark-brown background (vs. bright-orange light markings 
and blackish background) (Table 4). 

Phylogenetically and morphologically, Tylototriton kachinorum 
sp. nov. is most closely related to other species of Tylototriton 
inhabiting Myanmar (e.g., 7. shanorum and T. ngarsuensis) 
from the Shan Plateau and T. himalayanus from Nepalese 
Himalaya (Table 4). The new species can be readily distinguished 
from T. shanorum by having longer head in males (RHL 27.6 
vs. 22.4), greater internarial distance in both sexes (RIND 7.7— 
8.0 vs. 6.0—6.7), notably longer tail in both sexes (RTAL 120.5 
vs. 111.2 in males, 100.5 – 114.0 vs. 97.0—97.8 in females), 
wider vomerine tooth series in both sexes (RVTW 8.7—9.0 vs. 
6.1), longer vomerine tooth series in males (RVTL 12.8 vs. 
10.7) (Table 4), supratemporal bony ridges beginning at 
anterior corner of orbit (vs. supratemporal bony ridges beginning 
at loreal region), and wide, non-segmented vertebral ridge (vs. 
narrow, weakly segmented vertebral ridge). Tylototriton 
kachinorum sp. nov. has generally lighter coloration than T. 
shanorum: brownish ground color with orange-brown light 
markings (vs. more contrasting dark brown to black background 
color with yellow to reddish-brown light markings) and light 
markings present only on ventral tail ridge (vs. lateral sides of 
tail with light markings) (Table 4). 

Tylototriton kachinorum sp. nov. can be easily distinguished 
from T. ngarsuensis by having generally smaller body size in 
males (SVL 68.6+2.9 mm vs. 74.9-76.4 mm) and females 
(SVL 72.5-84.8 mm vs. 102.3 mm), comparatively longer 
head in males (RHL 27.6 vs. 24.0—-26.0) and females (RHL 
23.9—24.9 vs. 21.5), longer tail in males (RTAL 120.5 vs. 98.0- 
103.5), generally longer tail in females (RTAL 100.5—114.0 vs. 
104.6), and comparatively shorter vent length (in males RVL 
10.2 vs. 10.7—12.3, in females 5.6—6.8 vs. 8.0) (Table 4). The 
new species can be further distinguished from T. ngarsuensis 
by having truncate snout (vs. rounded snout), supratemporal 
ridges starting at anterior corner of orbit (vs. supratemporal 
bony ridges beginning posterior to orbit), non-segmented 
vertebral ridge (vs. weakly segmented), and 13-14 weakly 
distinct rib nodules (vs. 15 well-distinct rib nodules) (Table 4). 
Tylototriton kachinorum sp. nov. also has much lighter and 
duller colorationthan T. ngarsuensis: background color brown 
to light brown (vs. nearly black) with light orange-brown 
markings on rib nodules, parotids, and whole limbs (vs. no 
light markings on rib nodules or parotids, dark-yellow coloration 
present only on palms and soles) (Table 4). 

Morphologically, Tylototriton kachinorum sp. nov. most 
resembles its sister species, T. himalayanus from Nepal; 
however, it can be readily distinguished by the following 
morphological attributes: comparatively longer head in males 
(RHL 27.6 vs. 24.5), notably wider head in both sexes (RHW 
25.1—25.4 vs. 23.0- 23.6), shorter internarial distance in both 
sexes (RIND 7.7—8.0 vs. 8.2-8.4), and longer tail in both 
sexes (RTAL 100.5—120.5 vs. 98.0—100.6) (Table 4). The new 


species has generally smoother skin than T. himalayanus and 
can be distinguished from the latter species by having weakly 
distinct 13— 14 rib nodules (vs. large and prominent 16 rib 
nodules) and absence of lateral grooves on tail (vs. very 
distinct). Coloration of the new species is similar to T. 
himalayanus, but light markings are also present on labial 
regions, parotoids, rib nodules, and whole limbs (vs. no light 
markings on head and rib nodules, on limbs only on ventral 
surfaces) (Table 4). 


Etymology: The specific name “kachinorum” is a Latin 
adjective in the genitive plural (masculine gender), derived 
from the name of the Kachin people who inhabit the montane 
areas of northern Myanmar and adjacent territories (Kachin 
Hills), including the type locality of the new species. 


Recommended vernacular name: We recommend the 
following name in English: Kachin Crocodile Newt. 
Recommended vernacular name in Burmese (Myanmar) 
language: Kachin Yae Poke Thin. 


Conservation status: Tylototriton kachinorum sp. nov. is, to 
date, known from a single locality in the southern part of the 
Kachin Hills of northern Myanmar; the actual range of the new 
species is unknown. The new species is anticipated to inhabit 
elevations above 900 m a.s.l. on mountains surrounding the 
Indawgyi Lake valley and may be found in other parts of the 
Kachin Hills. Further research is required to estimate the 
actual distribution, population trends, and possible threats to 
the new species. Tylototriton kachinorum sp. nov. appears to 
be associated with montane forests and may be affected by 
growing anthropogenic pressure and forest destruction 
observed in different areas of Kachin State in Myanmar. Given 
the available information, we suggest Tylototriton kachinorum 
sp. nov. to be tentatively considered as a Vulnerable (VU) 
species following IUCN’s Red List categories (IUCN, 2001). 


DISCUSSION 


Our phylogenetic data largely confirmed the phylogeny of 
Tylototriton as presented in the earlier studies of Phimmachak 
et al. (2015), Khatiwada et al. (2015), Wang et al. (2018), and 
Grismer et al. (2018a). In particular, our data support the 
subdivision of Tylototriton into two major groups, traditionally 
regarded as subgenera: Tylototriton s. str. (clades 1 and 2) 
and Yaotriton (clades 3—5). Two other subgenera recently 
proposed by Fei et al. (2012), namely Qiantriton and 
Liangshantriton are nested within the radiation of Tylototriton 
s. str. (Figure 2) and possibly should not be regarded as 
independent subgenera or even genera (see Gong et al., 
2018; Hernandez, 2016). Our mtDNA-based genealogy is 
insufficiently resolved in a number of deeper tree nodes; 
however, it unambiguously suggests that clade 1 of Tylototriton 
S. str. is subdivided into three main subclades with prominent 
geographic structuring. The northernmost member of this 
group — T. kweichowensis, occurring in the Guizhou Plateau of 
China (Figure 1, samples 7—10), is distant from other members 
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of clade 1 and its phylogenetic position is unresolved. Most 
other species of clade 1, which occur in northern Indochina 
and Yunnan Province of China, form a monophyly (including T. 
verrucosus, T. shanjing, T. uyenoi, T. anguliceps, T. podichthys, 
T. pulcherrimus, T. yangi, and Tylototriton sp. 1) (Figure 1, 
samples 17—18 and 24—48). Finally, four species from the 
northern part of Myanmar and Himalaya form a distinct 
monophyletic group (including T. shanorum, T. ngarsuensis, T. 
himalayanus, and Tylototriton kachinorum sp. nov.). 

Studies on the taxonomic status of Tylototriton species in 
Myanmar have been delayed, even though the genus has 
been reported from the country for long time (see Gyi, 1969). 
Nishikawa et al. (2014) assigned Tylototriton from the Shan 
Plateau of eastern Myanmar to a new species, T. shanorum, 
whereas Grismer et al. (2018a) recently demonstrated the 
presence of divergent lineages of Tylototriton within different 
parts of the Shan Plateau and described a second species for 
the country — T. ngarsuensis. Contrary to the tree of Grismer 
et al. (2018a), our analysis suggests that T. shanorum is 
paraphyletic with respect to T. ngarsuensis and is subdivided 
into two lineages (Figure 2). However, it is worth noting that 
the second lineage of T. shanorum is based on data from a 
single specimen, reported by Nishikawa et al. (2014) 
(KUHE42348), which was obtained through the pet trade and 
presumed to come from Myanmar. Its assignment to T. 
shanorum is thus tentative and this specimen possibly 
represents a new lineage of Tylototriton from the Shan 
Plateau, yet undiscovered in the wild. Hence, considering the 
profound morphological differences between T. ngarsuensis 
апа Т. shanorum reported by Grismer et al. (2018a), though 
genetic divergence between these species was low (P=1.8% 
for ND2 gene), we consider that T. ngarsuensis indeed 
represents a distinct species of Tylototriton. 

Our description of Tylototriton kachinorum sp. nov. thus 
represents the third species of Tylototriton endemic to 
Myanmar. Tylototriton kachinorum sp. nov. is more closely 
related to T. himalayanus than to T. shanorum and T. 
ngarsuensis, occurring on the Shan Plateau in the eastern 
part of the country. This fact can be explained from a 
biogeographic viewpoint: the Kachin Hills represent the 
southernmost outcrop of the Great Himalaya ridge (Hadden, 
2008) and are isolated from the Shan Plateau by the 
Irrawaddy (Ayeyarwaddy) River valley, which may serve as an 
important biogeographic border for forest-dwelling taxa. 
Tylototritonhimalayanus, which can be found in Nepal and 
possibly along the Great Himalaya Ridge, was previously 
reported for northern Myanmar (Hernandez, 2016, 2018), but 
without any information on voucher specimens or reasons for such 
identification. There is a possibility that these records are based 
on the misidentification of Tylototriton kachinorum sp. nov. 

The present work indicates that Tylototriton diversity in 
Myanmar is still underestimated. Phimmachak et al. (2015) 
and Grismer et al. (2018a) reported sequences of Tylototriton 
cf. verrucosus from two regions of northern Myanmar: Sagaing 
Region (sample 18, see Figure 1) and the eastern part of 
Kachin State (sample 17, see Figure 1). Surprisingly, these 
populations were found to be distantly related to Tylototriton 
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kachinorum sp. nov., despite the geographical proximity to 
the southern part of Kachin Hills where the new species 
occurs. Sagaing and eastern Kachin Tylototriton populations 
(indicated herein as Tylototriton sp. 1) were reconstructed as 
members of the T. verrucosus species complex; however, they 
formed an mtDNA lineage clearly distinct from all currently 
recognized species. These data suggest that Tylototriton cf. 
verrucosus from northern Myanmar might represent a 
currently undescribed species and further morphological and 
phylogenetic studies are required to evaluate its taxonomic 
status. In agreement with the results of Wang et al. (2018), 
our analysis also shows deep subdivision and the presence of 
several lineages of possibly full-species status within T. 
asperrimus, T. dabienicus, and T. wenxianensis of the subgenus 
Yaotriton. 

Our study indicates that knowledge on amphibian diversity 
in montane regions of northern Myanmar is still far from 
complete and further diversity is likely to be revealed with 
additional survey efforts. A number of recent studies have 
shown that species diversity of amphibians and reptiles is 
widely underestimated across Myanmar (Grismer et al., 
2017a, 2017b, 2018a, 2018b, 2018c; Mulcahy et al., 2018), 
and it is likely that further surveys on this geographically 
complex and insufficiently studied region will lead to more 
discoveries. Recent economic development of Myanmar has 
led to increasing habitat loss and modification across the 
country (Li & Quan, 2017). Further intensified survey efforts 
and biodiversity assessments are urgently required for effective 
conservation management of yet unrealized herpetofaunal 
diversity in Myanmar. 
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ABSTRACT 


A new species of the genus Liurana Dubois, 1986 is 
described from Medog County, Tibet, China, based on 
morphological and molecular data. The new species 
can be differentiated from all other congeners by the 
following combination of characters: (1) head wider 
than long; (2) tympanum distinct and large; (3) hindlimb 
long, tibiotarsal articulation beyond tip of snout when 
adpressed; (4) belly with flat tubercles, cloacal region 
with small tubercles; (5) transverse bands distinctly on 
dorsal limbs, four bands on thigh and three on tibia; 
and, (6) dark brown marbled patterns or speckles on 
white belly. Here, we also discuss the distribution pattern 
of Liurana in the East Himalaya region, the role of the 
Yarlung Tsangpo River in the speciation and genetic 
isolation of congeners, the direct developmental mode of 


reproduction, and the two different ecotypes of the genus. 


Lastly, we provide conservation recommendations for the 
genus in southeastern Tibet. 


Keywords: Advertisement call; Biogeography; Ecology; 
Natural history; Tibet; Taxonomy 


INTRODUCTION 


As a wide-spread amphibian family in Southeast Asia, species 
of the family Ceratobatrachidae are distributed from the southern 
foothills of the Himalaya to the tropics of Southeast Asia (Yan et 
al., 2016). Members of the family are characterized as tropical 


Science Press 


specialists, with many species undergoing direct development 
without reliance on standing water bodies for breeding (Brown 
& Alcala, 1982; Brown et al., 2015). Within this family, frogs 
of the genus Liurana represent an understudied yet fascinating 
endemic group from the East Himalaya region. 

First considered as a subgenus of /ngerana by Dubois 
(1986), Liurana was established based on the type species 
Cornufer xizangensis Hu, 1987, with the subgenus later 
elevated to full genus based on morphological evidence (Fei 
et al, 1997). This taxonomic elevation is supported by 
recent phylogenetic studies, where Liurana was recovered as 
a distinct monophyletic clade from /ngerana, Platymantis, and 
Cornufer (Yan et al., 2016). However, despite research efforts 
on the higher-level systematic relationships of the genus Liurana, 
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little attention has been given to the species level taxonomy of the 
group in the Himalaya. To date, only a few studies have focused 
on species level taxonomy of the genus, with just a single study 
conducted in the last decade (Borah et al., 2013; Fei et al., 1997; 
Huang & Ye, 1997; Sichuan Biological Research Institute, 1977). 
As a result, all recognized species of Liurana are known from 
only a few vouchered specimens, and little is understood about 
their morphological variation and population structure. Based on 
limited studies, three species have been recognized in the genus 
to date, including L. alpina Huang, Ye, 1997, L. medogensis Fei, 
Ye, Huang, 1997, and L. xizangensis (Hu, 1987). 

During herpetological surveys in southeastern Tibet from 
2012 to 2016, 20 specimens of the genus Liurana were 
collected from Bomé and Medog Counties of the Nyingchi 
Prefecture, Tibet, China. Combining phylogenetic and 
morphological datasets, we describe here a new species of the 
genus from the tropical rainforest of Medog County, Nyingchi 
Prefecture, Tibet, China. Furthermore, we comment on the 
evolution, ecology, natural history, and conservation of the 
genus Liurana in China. 


MATERIALS AND METHODS 


Taxon sampling 

A total of 20 individuals of the recognized species of the genus 
Liurana were collected from different localities in the Bomé 
and Medog Counties, Nyingchi Prefecture, southeastern Tibet, 
China, and were comprised of 18 individuals of L. alpina, L. 
medogensis, and L. xizangensis and two individuals (one adult 
male and one adult female) of the new species (Figure 1; 
Appendix l). Following euthanasia, tissue samples were taken 


India 


Figure 1 Distribution of Liurana species in the East Himalaya 





and preserved in 95% ethanol, with the specimens then fixed in 
10% buffered formalin solution and transferred to 75% ethanol 
after fieldwork. All specimens were deposited in the Museum 
of the Kunming Institute of Zoology, Chinese Academy of 
Sciences (KIZ) (Appendix 11). 


Morphological data 

We confirmed the sex of each specimen by anatomical 
observation, with an incision made on the left side. All 
measurements were carried out by Ke Jiang using a digital 
caliper to the nearest 0.1 mm. Morphological characters and 
their measurement followed Fei et al. (2009) and included: 
snout-vent length (SVL); head length (HL), measured from 
posterior corner of mandible to tip of snout; head width (HW), 
measured at the angle of the jaw; snout length (SL), measured 
from tip of snout to anterior corner of eye; internarial distance 
(IND); interorbital distance (IOD), measured at the shortest 
distance between upper eyelid; maximum width of upper 
eyelid (UEW); eye diameter (ED), measured as the horizontal 
diameter of eye); tympanum diameter (TD), measured as the 
maximum horizontal diameter of tympanum); length of lower 
arm and hand (LAHL), measured from the elbow joint to the tip 
of the longest finger; largest diameter of lower arm (LAD); hand 
length (HAL), measured from the base of the outer metacarpal 
tubercle to the tip of finger III; femur length (FML), measured as 
the linear distance between the insertion of the leg to the knee 
joint; tibia length (TL), measured as the linear distance between 
the knee joint and tibiotarsal articulation; length of tarsus and 
foot (TFL), measured from the tibiotarsal articulation to the tip 
of toe IV; and foot length (FL), measured from the base of the 
inner metatarsal tubercle to the tip of toe IV. 
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blue localities represent reported records of Liurana cf. medogensis from Southern Tibet (see discussion for details). 
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In addition to the newly obtained data, morphological data 
of congeners were also obtained from published literature for 
comparison (Borah et al., 2013; Fei et al., 1997, 2009; Huang 
& Ye, 1997). 


Molecular analysis 
Genomic DNA was extracted from tissue samples using 


standard phenol-chloroform protocols (Sambrook et al., 1989). 


Fragments of a single mitochondrial DNA locus (cytochrome 
c oxidase subunit |, CO/) and three nuclear loci, including 
recombination activating protein 1 (Rag1), tyrosinase (Tyr), 
and rhodopsin (Rhod), were targeted and amplified following 
published primers and protocols (Che et al., 2011; Yan et 
al., 2016). The products were purified with a Gel Extraction 
Mini Kit (Watson Biotechnologies, Shanghai, China) and 
sequenced on an ABI 3730xl DNA automated sequencer 
(Applied Biosystems, Uk). 

Additional sequences of congeners and closely related 


outgroups were obtained through GenBank (Appendix Il). 


Sequences were aligned using MUSCLE implemented in 
Geneious R10.0.6. Uncorrected genetic distances of the 


mitochondrial locus COI were calculated using PAUP* v4.0b10. 


To assess phylogenetic congruence between the mitochondrial 
and nuclear data, the phylogeny for each gene was inferred 
independently using Bayesian analyses. As по strongly 
supported incongruences between mitochondrial and nuclear 
data were observed visually, two datasets were concatenated 
for the final analyses. 

Phylogenetic relationships were inferred using both 
partitioned Bayesian (MrBayes v3.2.1; Ronquist & 
Huelsenbeck, 2003) and maximum likelihood analyses 
(RAXML VI-HPC v8.2.10; Stamatakis, 2014). For Bayesian 
analyses, all protein-coding genes were partitioned by codon 
positions, and the best models of nucleotide substitution were 
selected for each partition by the Akaike Information Criterion 
(AIC), as implemented in jModelTest v2.1.10 (Darriba et al., 
2012; Guindon & Gascuel, 2003). A rate multiplier model 


was used to allow substitution rates to vary among subsets. 


Default settings were used for all other model parameters. 
Two independent Markov chain Monte Carlo analyses were 
run, each with four Metropolis-coupled chains, a melting 
temperature of 0.02, and an exponential distribution with a rate 
parameter of 25 as the prior branch lengths (Marshall, 2010). 
All Bayesian analyses were run for 6 000 000 generations, with 
parameters and topologies sampled every 3000 generations. 
Stationarity and convergence were assessed with TRACER 
v1.6.0 (Rambaut et al., 2013). 

Partitioned maximum likelihood analyses were conducted on 
the concatenated dataset using the same partitioning strategy 
as for Bayesian analysis. The more complex model (GTR«I) 
was applied for all subsets (Table 1), with 1000 replicate ML 
inferences. Each inference was initiated with a random starting 
tree, and nodal support was assessed with 1000 bootstrap 
pseudoreplicates (Stamatakis et al., 2008). 

The resulting phylogenetic trees were rooted using the clade 
containing /ngerana and Limnonectes, following recent studies 


on phylogenetic relationships of the focal group (Yan et al., 2016). 


Table 1 Partition strategies and best evolutionary models 
selected for each partition 


Gene Codon Fragment length (bp) Model selected 








col 1st 190 GTR 
2nd 190 НКҮ+Г 
зга 190 НКҮ 
Rag1 1st 184 F81 
2nd 184 HKY 
3rd 184 HKY 
Rod 1st 105 JC 
2nd 105 F81 
3rd 105 GTR 
Tyr 1st 177 ОТН+Г 
2nd 177 HKY«T 
3rd 177 HKY«T 
RESULTS 
Morphology 


Morphometric variation of the examined Liurana species is 
summarized in Table 2. The unidentified specimens of Liurana 
were morphologically most similar to L. medogensis and could 
be differentiated from all recognized species by having a 
relatively wider head (HW/HL >100% vs. <100%) and longer 
hindlimbs, with tibiotarsal articulation reaching beyond tip of 
snout when adpressed (vs. reaching anterior corner of eye 
only). Furthermore, the unidentified Liurana specimens from 
Medog possessed distinct tubercles near the cloaca, which are 
absent or indistinct in all recognized congeners except for L. 
xizangensis. 

Coloration and ornamentation were highly variable among 
the examined specimens of L. alpina and L. xizangensis, 
ranging from uniform bright orange reddish to brownish gray 
with blackish speckles (Figure 3). For L. medogensis and 
the two unidentified individuals from Medog, the coloration 
and ornamentation were less variable. A single individual 
of L. medogensis (KIZ05587) possessed a broad, yellowish 
dorsal vertebral stripe from snout to vent, whereas the other 
individuals of the same species exhibited a much darker 
vertebral stripe in light reddish brown. 


Phylogeny 

The phylogenetic placement of the genus Liurana is similar to 
previous results (Yan et al., 2016), with the genus recovered as 
monophyletic with strong support (Bayesian posterior probability 
(PP)=1.00; ML bootstrap support (BS)=100) (Figure 2). Within 
the genus, L. xizangensis and L. alpina form a monophyletic 
group with strong supports (PP=0.99, BS=100), which is sister 
to L. medogensis (PP=1.00, BS=100). The two unidentified 
individuals (Liurana sp.) collected from Medog County form a 
distinct, monophyletic clade (PP=1.00, BS=100), which is basal 
with respect to all other Liurana congeners (Figure 2). 
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For the obtained 570 bp fragment of CO/ data, the 
uncorrected genetic distance between the two unidentified 
Liurana individuals is 0.88% only, while such distance 
is 22.19%-22.22% to its morphologically most similar 
congener, L. medogensis, 21.67%—22.63% to L. alpina, and 
20.19%—21.30% to L. xizangensis (Table 3). The observed 
intraspecific variation among recognized Liurana species is 
0%-0.183% within conspecifics of L. alpina. For L. xizangensis, 
the greatest genetic divergence (3.12%) is observed between 
individuals from Tongmai, Bomé County, Tibet, China (Locality 
#1) and 62K, Medog County, Tibet, China (Locality #2), which 
are on different sides of the Yarlung Tsangpo River (Figure 1). 

Because the two unidentified specimens of Liurana 


from Medog possess a suite of distinct morphological 
characteristics, high genetic divergence from congeners, and 
monophyletic and distinct phylogenetic position, we here 
describe them as a new species. 


Taxonomy 


Liurana vallecula sp. nov. Jiang, Wang, Wang, Li, and Che 
Figures 1—5. 


Holotype: КІ2014083, adult female, from the side of Yarlung 
Zangbo River near Xirang (N29.1776°, E95.0057°, elevation 
550 m a.s.l.; Figure 3), Beibeng, Medog County, southeastern 
Tibet, China. Collected by Cheng Li on 17 April 2016. 


Table 2 Morphological comparisons for four species of the genus Liurana 








Liurana vallecula sp. nov.  L. alpina L. medogensis L. xizangensis 
SVL (mm) 14.6 (19) 23.2-24.9 (33) 13.1-19.0 (37) 20.6-22.5 (97) 
20.4 (19) 25.5 (19) 29.4-30.5 (22) 
HL/HW <100% >100% >100% >100% 
Tibiotarsal articulation Beyond tip of snout Reaching anterior corner of eye Reaching anterior corner of eye Reaching anterior corner of eye 
Flat tubercles on belly Absent Absent Absent Present 
Tubercles around cloaca Present Absent Indistinct Present 


Ventral pattern Thin, marbled-patterns 


Morphological abbreviations are listed in the methods. 


Thin, marbled-patterns 


Thick, vermiculated stripes Thin, marbled-patterns 


Table 3 Uncorrected genetic distances of the CO/ fragment obtained (570 bp) among Liurana species 





1 L. sp. nov. YPX47504 


2 L. sp. nov. YPX47527 0.00877 – 

3 L. medogensis KIZ 010955 0.22191 0.222 = 

4 L. alpina KIZ 07357 0.22029 0.21666 0.13394 — 

5 L. alpina KIZ 011140 0.22632 0.22281 0.13207 0.00178 
6 L. alpina KIZ 011141 0.22216 0.21853 0.13211 0.00183 
7 L. xizangensis KIZ 06707 (Locality #2) 0.2129 0.21297 0.14312 0.09541 
8 L. xizangensis KIZ 09956 (Locality #2) 0.2129 0.21297 0.14312 0.09541 
9 L. xizangensis KIZ 012705 (Locality #2) 0.2129 0.21297 0.14312 0.09541 
10 L. xizangensis KIZ 012706 (Locality #2) 0.2129 0.21297 0.14312 0.09541 
11 L. xizangensis KIZ 011104 (Locality #2) 0.2129 0.21297 0.14312 0.09541 
12 L.xizangensis KIZ 012707 (Locality #2) 0.21107 0.21114 0.14495 0.09725 
13 L. xizangensis KIZ 012704 (Topotype) 0.20372 0.20379 0.12661 0.08624 
14 L. xizangensis KIZ 014046 (Locality #1) 0.20187 0.20195 0.12844 0.08807 


0 
0.09345 
0.09345 
0.09345 
0.09345 
0.09345 
0.09535 
0.08436 
0.0862 


0.09358 – 

0.09358 0 - 

0.09358 0 0 - 

0.09358 0 0 0 = 

0.09358 0 0 0 0 = 

0.09541 0.00183 0.00183 0.00183 0.00183 0.00183 – 

0.0844 0.02752 0.02752 0.02752 0.02752 0.02752 0.02936 - 
0.08624 0.02936 0.02936 0.02936 0.02936 0.02936 0.03119 0.00183 – 





Paratype: КІ2014106, adult male, from the side of Baimaxilu 
River (a tributary of the Yarlung Zangbo River) near Maniweng 
(N29.2667°, E95.1667°, elevation 850 m a.s.l.), Beibeng, 
Medog County, southeastern Tibet, China. Collected by Ke 
Jiang and Yu-Fan Wang on 29 April 2016. 


Diagnosis: Liurana vallecula sp. nov. is assigned to the 
genus Liurana by its phylogenetic position and the following 
morphological characters: (1) body size small (SVL 14.6—20.4 mm, 
n=2); (2) tips of fingers and toes not expanded, (3) grooves absent 
on tips of fingers and toes; (4) webbing absent on all digits; (5) 
metacarpal tubercles and metatarsal tubercles absent; (6) tarsal 
fold absent; and (7) vocal sac and vocal sac openings absent. 
Liurana vallecula sp. nov. can be distinguished from 
all currently recognized congeners by a combination of the 
following morphological characters: (1) head wider than long 
HW/HL 108%—113%; (2) tympanum distinct, large; (3) hind limbs 
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long, with tibiotarsal articulation reaching beyond tip of snout 
when adpressed; (4) flat tubercles present on ventral surface of 
body; (5) small tubercles present on cloacal region; and (6) dark 
brown marbled pattern present on ventral surface of body. 


Description of holotype: Body size small, SVL 20.4 mm. 
Head wider than long (HW/HL=108%); snout rounded, slightly 
projecting over lower jaw; canthus rostralis distinct, relatively 
sharp; nostril closer to tip of snout than to eye; loreal region 
slightly concave and oblique; eye relatively large (ED/HL=40%); 
internarial distance larger than interorbital distance and width 
of upper eyelid; tympanum oval, concealed under skin, 
with distinct, elevated rim, slightly less than half of eye 
diameter (TD/ED=45%); supratympanic fold distinct, extending 
posterior-inferiorly to area above forelimb insertion; pair of 
distinct skin folds present on head from posterior orbit dorsally 
to pectoral region, symmetrical along vertebral line, extending 
toward medial line posteriorly. 


Ingerana tenasserimensis 
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Figure 2 Phylogenetic relationships among Liurana species based on maximum likelihood and Bayesian analyses of one 


mitochondrial gene (СОГ) and three nuclear genes (Rag1, Rod, and Tyr) 
Maximum likelihood bootstrap and Bayesian posterior probability values are given at all nodes (in such order respectively), except the internal nodes within L. 


xizangensis from Locality 42, which have short branch lengths and resemble polytomy. 





Figure 3 Holotype of Liurana vallecula sp. nov. in life (adult female, KIZ014083) 
A: Dorsolateral view; B: Ventral view; C: Ventral close-up of hand; D: Ventral close-up of feet. Photos by Yu-Fan Wang. 
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Forelimb slender; forearm and hand length less than half of 
body length (LAHD/SVL=44%); fingers compressed vertically, 
tips rounded but not expanded, transverse grooves absent, 
relative finger lengths I«II«IV«III; subarticular tubercles absent; 
three metacarpal tubercles present, flat and indistinct. Hind 
limb slender, tibiotarsal articulation reaching beyond tip of 
snout when adpressed; heels overlapped when flexed and 
held perpendicular to body; tibia length larger than half 
of body length (TL/SVL=58%); toes compressed vertically, 
tips slightly expanded, transverse grooves absent; relative 
toe lengths I<ll<V<lll<IlV; toe webbing absent; subarticular 
tubercles indistinct; inner metatarsal tubercle oval, indistinct; 
outer metatarsal tubercle indistinct; tarsal fold absent. 

Relatively weak, discontinuous folds present on dorsolateral 
body from shoulder to about two thirds of trunk on each 
side of body; single skin fold present along vertebrate from 
snout to vent, much weaker than discontinuous folds laterally; 
dorsal surface relatively rough, tubercles randomly scattered 
on dorsal and lateral head, body, and limbs as well as around 
cloaca; tubercles much finer on ventrolateral region of body. 
Ventral head, body, and limbs mostly smooth, except several 
flat tubercles on base of ventral thigh and small tubercles 
on ventral surface of tarsus and metatarsus. Vomerine 
teeth absent; tongue large, elongated oval, deeply notched 
posteriorly, small papillae scattered on tongue. 


Coloration of holotype in life: The dorsal surfaces of the 
head and body are reddish brown. Dark brown streaks and 
marble patterns are present on dorsal head, body, and limbs, 
including a transverse streak between orbit on dorsal surface of 
head, a X-shaped pattern on pectoral region of dorsum, irregular 
marble patterns on lateral body, and transverse streaks across 
dorsal limbs (more distinct on the dorsal hind limbs). The lower 
parts of canthus rostralis and supratympanic fold are blackish 
brown. The ventrolateral surface of the body is dark brown, with 
small white spots scattered across. Ventral surfaces of limbs, 
head, and body are light grayish brown; white marble patterns are 
present on ventral surfaces of head, body, and ventral forelimbs, 
and ventrolateral surfaces of thigh. The white marble patterns are 
finer and much smaller on ventral head comparing to ventral body, 
giving a speckle impression; white marble patterns on ventral 
body are mostly interconnected. A few smaller white spots are 
also present on ventral thigh, tibia, and femur. 


Coloration of holotype in preservative: Ornamentation 
patterns remain after preservation. | However, coloration 
changes after preservation, include: (1) dorsal surfaces of 
head, limbs, and body become grayish brown, with dark gray 
patches; (2) lower parts of canthus rostralis and supratympanic 
fold, transverse bands on dorsal limbs, as well as lower part 
of lateral body become dark gray; (3) ventral surfaces of throat 
and limbs become brown, with grayish white spots. 


Variation: Morphometric variations of the type series are 
shown in Table 4. Most external morphological characters 
are identical between the two individuals, but the paratype is 
smaller than the holotype (SVL 14.6 mm in paratype male vs. 
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20.4 mm in holotye female), as well as in having rather paler dorsal 
coloration and more dense spots on ventral body. No secondary 
sexual characters, such as vocal sac or nuptial pad, are observed 
in the paratype male, but a single black testicle was observed 
on the left side of the male paratype, which is oval shaped and 
relatively large, with a longitudinal length of about 1.5 mm. 


Ecological and natural history notes: Liurana vallecula 
sp. nov. is a terrestrial, leaf-litter specialist, inhabiting the 
forest floor of tropical broad-leaf forest at low elevations (below 
1000 m a.s.l. near Yarlung Zangbo River and its immediate 
tributaries. The female holotype had about five immature eggs 
in the left ovary, which were well developed and relatively large. 


Distribution: Currently the new species is known only from the 
type localities of Xirang and Maniweng of Beibeng, Medog County, 
Nyingchi Prefecture, Tibet, China. The new species likely inhabits 
other nearby regions in southern Tibet (see Discussion below). 


Etymology: The specific epithet of the new species, "vallecula" 
means “valley inhabitor’, in reference to the habitat of this 
species in the lower river valley of Yarlung Zangbo Grand 
Canyon. We suggest Valley Papilla-tongued Frog as its English 
common name and He Gu She Tu Wa (in| 4 2 E) as its 
Chinese common name. 





Comparisons: Liurana vallecula sp. nov. differs from 
the three congeners by having a much wider head (HW/HL 
>100% vs. <100%) and longer hindlimbs (tibiotarsal articulation 
reaching beyond tip of snout when adpressed vs. reaching only 
anterior corner of eye when adpressed). In addition, Liurana 
vallecula sp. nov. can be further distinguished from L. alpina 
and L. xizangensis by having a smaller body size (SVL 14.6 
mm in male, 20.4 mm in female vs. SVL 23.2—24.9 mm іп 
males, 25.5 mm in female for L. alpina; 20.6—24.5 mm іп 
males, 29.4—30.5 mm in females for L. xizangensis); and from 
L. medogensis by its different ventral pattern (small, marbled 
patterns or speckles vs. broad, vermiculated stripes). 


DISCUSSION 


Liurana medogensis from Southern Tibet 

Recently, Borah et al. (2013) reported L. medogensis from 
the eastern part of Southern Tibet, which they recognized as 
Limnonectes (Taylorana) medogensis in the study (localities 
1 and 2 in light blue; Figure 1) (all abbreviations of “L.” in 
this paragraph refer to Liurana, not Limnonectes). Although 
one of the reported specimens (ZSI a11549) resembles the 
external appearance of L. medogensis based on the figure 
in the manuscript (Borah et al., 2013; Figure 1), according to 
the description, it possesses morphological characteristics that 
differ from the diagnosis of the genus Liurana, including having 
a distinct vomerine ridge, prominent vomerine teeth, and 
rudimentary webbing on toes (vs. vomerine ridge, vomerine 
teeth, and toe webs absent in Liurana) (Fei et al., 2009, 2010; 
Yan et al., 2016). Furthermore, this specimen was collected at 
a much higher elevation (2 000—2 500 m a.s.l.) than the known 


range for L. medogensis (about 1 400 m a.s.l.). For the second 
reported specimen of L. cf. medogensis (BMHE a0081), based 
on the description provided by the authors, it matches the 
morphological diagnosis of our new species, Liurana vallecula, 
including having a wide head (HW/HL>100%), but lacks the 
broad, vermiculate patterns on the ventral surface of the 
body (Borah et al., 2013). As we cannot examine these 
specimens or obtain their genetic data, we cannot confirm the 
taxonomic status of these two specimens reported by Borah 
et al. (2013). Future taxonomic studies are needed to gain a 
better understanding of Liurana diversity in this region. 


Table 4 Measurements (in mm) of type series of Liurana 
vallecula sp. nov. 


























Number Status Sex SVL HL HW 51 IND IOD ЧЕМ 

2014083 Holotype Female 24 73 79 3.2 30 22 2.0 

Rati 

айоо - 372 387 156 145 107 98 

SVL (%) 

Z014106 Paratype Male 146 56 63 25 23 1.5 155 
Ratio to SVL (%) - 38.1 43.2 17.5 15.7 10.4 10.0 

Number ED TD LAHL LAD HAL FML TL TFL FL 

2014083 2.9 1.3 9.0 15 51 10.9 11.8 17.4 11.6 
Ratio to SVL (%) 14.1 6.4 44.1 73 247 533 576 853 57.0 

2014106 2.0 1.0 6.3 14 3.3 7.7 85 12.71 
Ratio to SVL (96) 13.7 43.6 7.4 22.4 52.7 58.4 76.8 49.0 43.6 





See methods section for abbreviations. —: Not applicable. 


Ecotypes 

Liurana diversity can be divided into two major ecotypes: 
the alpine ecotype that inhabits cool, moist mixed forests at 
2000-3000 m a.s.l. (including L. alpina and L. xizangensis), 
and the tropical ecotype that inhabits low-elevation tropical 
rainforest below 2000 m а.5.1. (including L. medogensis 
and Liurana vallecula). For the first ecotype, coloration of 
individuals is highly variable, ranging from uniform bright 
reddish orange to marbled purplish gray. Frogs of this ecotype 
live under thick layers of moss on fallen tree trunks or rocks 
along the forest edge. The second ecotype consists of leaf litter 
specialists that inhabit the forest floor under pristine tropical 
rainforest. Coloration of this ecotype is much more cryptic, 
ranging from a single, wide, brownish orange dorsal stripe to 
dark brown with darker marbling (Figures 4, 5). 

Even though L. medogensis and Liurana vallecula both 
belong to the second ecotype, they do not form a monophyletic 
group according to our genetic data (Figure 2). Future 
ecological studies are needed to further differentiate the 
ecological niche of each species and investigate the evolution 
of ecotypes in the genus Liurana. 


Reproductive biology 

Although the reproductive biology, particularly the direct 
developing reproductive mode, has long been documented 
within the family Ceratobatrachidae (Brown & Alcala, 1982; 
Siler et al., 2010), little is known about the reproductive biology 
of Liurana specifically. Considering the Tibetan Liurana species 
as Platymantis at that time, Hu et al. (1987) first commented on 
the reproductive mode of the Tibetan species and suggested 


that they may go through direct development, similar to other 
Platymantis species. It was only until 2010 when the eggs 
of Liurana species in the alpine ecotype from Tibet were first 
collected (Li et al., 2010). According to Li et al. (2010), eggs of 
the unidentified Liurana species from 62K in the Medog County 
(confirmed as L. xizangensis according to Yan et al. (2016) and 
the present study) were large (with a diameter of about 3.5 mm) 
and the clutch size was small. These egg characteristics are 
similar to those of direct development in the same family (Siler 
et al., 2010). Similar to Li et al. (2010), we also observed eggs 
in a female L. xizangensis (KIZ014153) from the same locality, 
where the female displayed 14 and 16 yellowish white eggs in 
the left and right oviduct, respectively, though not all were fully 
developed (Figure 6). 

More recently, Borah et al. (2013) provided information 
regarding the reproductive biology of the tropical ecotype 
of Liurana from Southern Tibet. Embryos displayed 
characteristics of the direct developing species in the family 
Ceratobatrachidae, including having a few embryos, large 
individual embryo size, and whitish coloration (Borah et al., 
2013). Based on these similarities, the authors claim that 
this Liurana species reproduces through direct development 
without the larval stage (Borah et al., 2013). However, such 
conclusions on direct development in Liurana are still based 
on indirect inferences of egg characteristics. Therefore, we 
recommend future studies focus on the reproductive biology of 
the genus to clarify the reproductive mode of Liurana. 


Acoustic signals 

In addition to the differences in general habitat preferences and 
external morphology, the two ecotypes of Liurana also differ in 
their acoustic signals. In fact, extensive acoustic signals have 
been observed in the alpine ecotype only (Fei et al., 2009; Hu et 
al., 1987; present study). Such differences in acoustic signals 
may be explained by the relative cost-benefit ratios of calling 
in each specific environment (Ryan, 1988). Future ecological 
and behavioral studies are needed to confirm this hypothesis 
in the field and investigate the communication strategies of the 
tropical ecotype. 


Evolution 

Frogs of the genus Liurana are endemic to eastern Himalaya, 
and all are found in the Yarlung Zangbo Drainage in Southern 
Tibet (Figure 1). Based on the available distribution data, L. 
xizangensis is distributed on the eastern side of the Yarlung 
Zangbo River, whereas L. alpina is found on the western 
side of the river. As the two species are found in cool 
environments at high elevation only (22000 m а.5.1.), the hot 
tropical valley of the Yarlung Zangbo River may serve as a 
dispersal barrier, thereby shaping the genetic diversity and 
facilitating the speciation processes of the two congeners. This 
hypothesis is partly supported by our genetic data, where 
closely distributed populations from two sides of the river 
possessed higher genetic divergence than further distributed 
populations from the same side (Table 3). Future studies 
should expand population sampling of the two species along 
both sides of Yarlung Zangbo River to examine its role in the 
evolution of Liurana. 
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L. alpina (C and D), L. 


Figure 5 Habitat of Liurana vallecula sp. nov. (A) (Xirang, elevation 550 m a.s.l.), L. medogensis (B) (Xigong Lake, elevation 1 300 m a.s.l.), 
L. alpina (C) (Dayandong, elevation 3 000 m a.s.l.), and L. xizangensis (D) (62K, elevation 2800 m a.s.l.) in Medog County, southeastern 
Tibet, China, respectively. Photos by Ke Jiang, Kai Wang, and Shuai Wang 
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Figure 6 Gravid female (A) and developing eggs through transparent abdomen (B) of Liurana xizangensis from 62K, Medog County, 


Tibet, China. Photos by Ke Jiang 


Conservation 

As micro-endemic habitat specialists, Liurana species are 
threatened by habitat destruction in southern Tibet. Based 
on our continuous field surveys since 2012, considerable 
habitat destruction has been observed at 62K in Medog 
County, which is one of only three known localities for L. 
xizangensis. Unregulated infrastructure developments have 
destroyed the mossy fields along the forest edges, streams, 
and wetlands, which constitute the core habitats not only for 
L. xizangensis, but also other micro-endemic anuran species 
such as Scutiger wuguanfui and S. spinosus (Jiang et al., 
2016). Similarly, continuous tourist development and road 
construction along the hiking trail of Medog pose serious 
threats to habitat at the only known locality of L. alpina. 
Therefore, we recommend that local authorities and regional 
governments take habitat conservation into account when 
making developmental decisions, and we urge law enforcement 
agencies to enforce the existing environmental regulations 
of construction projects in the region, particularly in Medog 
County. 
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APPENDIX | 
Specimens examined for comparison. Museum abbreviations follow those 
from the methods section. 


Liurana alpina (n=4): KIZ011140—42 (males), KIZ07358 (Dayandong (type 
locality)) Medog County, Tibet, China. 


Liurana medogensis (п=3): КІ2010955, 05886, 05587 (males), (Xigong Lake 
(type locality)) Medog County, Tibet, China. 


Liurana xizangensis (n=11): KIZ012704, 014046 (males), (Tongmai (close to 
type locality “Yi'ong”, =Yigong)) Bomé, Tibet, China; КІ206707, 09954—56, 
011105—06, 012706 (males); KIZ09953, 014153 (females), 62K, Medog 
County, Tibet, China. 


APPENDIX II GenBank sequences used and accession Nos. of novel sequences 


Gene and GenBank accession No. 























Genus Species Voucher or tissue number Locality col Rag1 Rod Tyr 
Alcalus tasanae CAS232349 Myanmar - KU243096 KU243106 U243116 
Alcalus tasanae CAS247243 Myanmar = KU243097 KU243107 U243117 
Alcalus tasanae THNHM20534 Myanmar - KU243098 KU243108 U243118 
Ceratobatrachus guentheri VUB1017(SR5543) Solomon Islands AY883979 DQ347272 00347391  DQ347179 
Cornufer boulengeri BPBM22329 Papua New Guinea НО844999 KP298309 – P298385 
Ingerana tenasserimensis | CAS205064/ TADP918 Myanmar/Thailand KR087736 00347258 AY322236 P298327 
Limnonectes limgorgi VUB1218 Laos = DQ347286 09347407 009347194 
Liurana alpina КІ2011140 Medog, Tibet, China МК462138 KU243094 KU243104 U243114 
Liurana alpina КІ2011141 Medog, Tibet, China МК462139 KU243095 KU243105 U243115 
Liurana alpina KIZ07357 Medog, Tibet, China | MK462137  — - _ 

Liurana medogensis KIZ010955 Medog, Tibet, China | MK462136  — KU243103 U243113 
Liurana vallecula KIZ014083 Medog, Tibet, China | MK462134 MK315115 МК462150 МК462155 
Liurana vallecula КІ2014106 Medog, Tibet, China МК462135  — MK462151 MK462156 
Liurana xizangensis KIZ06707 Medog, Tibet, China МК462140 KU243092  KU243101 U243111 
Liurana xizangensis КІ209956 Medog, Tibet, China МК462141 – - - 

Liurana xizangensis КІ2011104 Medog, Tibet, China МК462144  — - _ 

Liurana xizangensis KIZ012705 Medog, Tibet, China МК462142  — _ _ 

Liurana xizangensis KIZ012706 Medog, Tibet, China МК462143 — - _ 

Liurana xizangensis KIZ012707 Medog, Tibet, China | MK462145  — - _ 

Liurana xizangensis KIZ012704 Medog, Tibet, China | MK462146  — - _ 

Liurana xizangensis KIZ014046 Medog, Tibet, China | MK462147  — - _ 
Platymantis hazelae TNHC62160 & CMNH-RSK3918 Phillipines = DQ347248 DQ347369 DQ347153 

—: Not available. 
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ABSTRACT 


Insulin is a key hormone for the regulation of 
metabolism in vertebrates. Insulin is produced by 
pancreatic islet cells in response to elevated glucose 
levels and leads to the uptake of glucose by tissues 
such as liver and adipose tissue to store energy. 
Insulin also has additional functions in regulating 
development. Previous work has shown that the 
proglucagon gene, which encodes hormones counter 
regulating insulin, is duplicated in teleost fish, and that 
the peptide hormones encoded by these genes have 
diversified in function. | sought to determine whether 
similar processes have occurred to insulin genes in 
these species. Searches of fish genomes revealed an 
unexpected diversity of insulin genes. A triplication of 
the insulin gene occurred at the origin of teleost fish, 
however one of these three genes, insc, has been lost 
in most teleost fish lineages. The two other insulin 
genes, insa and insb, have been retained but show 
differing levels of selective constraint suggesting that 
they might have diversified in function. Intriguingly, 
a duplicate copy of the insa gene, which | named 
insab, is found in many fish. The coding sequence 
encoded by insab genes is under weak selective 
constraint, with its predicted protein sequences losing 
their potential to be processed into a two-peptide 
hormone. However, these sequences have retained 
perfectly conserved cystine residues, suggesting that 
they maintain insulin’s three-dimensional structure and 
therefore might modulate the processing and secretion 
of insulin produced by the other genes. 


Keywords: Insulin; Teleost fish; Gene duplication; 
Adaptive evolution; Gene loss 


Science Press 


INTRODUCTION 


Fish have been important contributors to our understanding 
of human biology and disease (Lieschke & Currie, 2007; 
MacRae & Peterson, 2015), especially in endocrinology 
(Conlon, 2000a). In vertebrates, insulin is a hormone produced 
by the beta-cells of pancreatic islets in response to increased 
levels of blood glucose, which induces the uptake of glucose 
by tissues such as liver and adipose tissue for storage, and 
thus is a key regulator of glucose metabolism (Róder et al., 
2016). Deficiencies in insulin production and/or signaling 
leads to diabetes in humans and other animals (Róder et al., 
2016; Weiss, 2009). Many aspects of insulin function have 
been learned from extensive studies in diverse species of fish 
(Caruso & Sheridan, 2011; Conlon, 2000b; Polakof et al., 2012). 
Fish have also been explored in the development of potential 
treatments for diabetes. Xenotransplantation of fish islets, using 
both wild-type and humanized insulin, has been considered due 
to the relative ease of isolation of islets from pancreatic tissue in 
fish (e.g., tilapia, Oreochromis niloticus) compared to mammalian 
sources (Wright et al., 2014). Insulin, however, in addition to 
its role in the regulation of metabolism, also has other roles 
in vertebrate biology, including some in development (Chan & 
Steiner, 2000; Hernández-Sánchez et al., 2006). 

Insulin is a polypeptide hormone composed of two peptides, 
A- and B-chains of about 20 and 30 amino acids, respectively, 
which are held together by disulfide bridges (Conlon, 2001; 
Steiner et al., 1985, 2009; Weiss, 2009). Insulin sequences 
have been characterized, mostly at the protein level, in a large 
number of vertebrate and fish species (Caruso & Sheridan, 
2011; Conlon, 2000b, 2001). Typically, vertebrate species have 
a single copy of the insulin gene in their genomes (Nishi & 
Nanjo, 2011; Steiner et al., 1985), however, multiple insulin 
genes have been characterized in some species, such as 
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rodents (Long et al., 2013; Soares et al., 1985; Wentworth et al., 
1986) and some fish (Caruso et al., 2008; Caruso & Sheridan, 
2011, 2014; Irwin, 2004). Within fish, it appears that only 
one of these insulin genes, which produces the hormone that 
regulates blood glucose levels, is expressed in the pancreas 
(Caruso & Sheridan, 2011; Irwin, 2004), with the exception of 
the multiple insulin genes produced by the very recent genome 
duplication in salmonid fish (Caruso et al., 2008; Caruso & 
Sheridan, 2014). While the origin of the two insulin genes 
found in some fish (e.g., zebrafish) was shown to be early 
in teleost fish evolution, whether this ins2 gene has been 
retained in the genomes of diverse fish is unclear (Caruso & 
Sheridan, 2011; Irwin, 2004). The fish ins2 gene has been 
poorly characterized. The initial identification of this gene in 
zebrafish (Danio rerio) and takifugu (Takifugu rubripes) only 
showed expression in the embryo of zebrafish (Irwin, 2004). A 
subsequence study of the expression of the two insulin genes 
in zebrafish suggested that ins2 potentially has a function as a 
growth and neurotrophic factor during development (Papasani 
et al., 2006). More recently, studies of the tilapia (O. niloticus) 
ins2 gene showed that it had negligible expression in the 
pancreas, thus likely is not a major contributor to the regulation 
of glucose metabolism and would not need to be silenced to 
allow the xenotransplantation of tilapia islets as a treatment for 
diabetes in humans (Hrytsenko et al., 2016). 

Changes in the function of insulin can be paralleled by 
changes in the function of glucagon, the hormone that counters 
the action of insulin (Seino et al., 1988). Recently it has been 
shown that duplicated proglucagon genes are widespread in 
teleost fish, and that the functions of the proglucagon-derived 
peptides encoded by these genes have diversified (Irwin & 
Mojsov, 2018). Here | have taken advantage of the large 
number of fish genomes that have been characterized in the 
past few years (Ravi & Venkatesh, 2018) to determine the 


distribution of insulin genes in the genomes of ray-finned fish. 


Analysis of these sequences allows us to begin to understand 
the origin of these genes, including the measurement of 
the selective forces acting upon them, which might provide 
evidence for differences in function. Surprisingly, | found a 
greater diversity in the numbers of insulin genes than expected, 
with the identification of a third insulin paralog, as well as 
lineage-specific duplicates of some of these paralogs. My 
analyses suggest that all three fish insulin paralogs encode 
functional protein products. 


MATERIALS AND METHODS 


Data collection 
Genome sequences databases maintained by the National 


Center for Biotechnology Information (NCBI: https://www.ncbi. 
nlm.nih.gov/genome/gdv/) and Ensembl (http://www.ensembl. 


org) were searched in March/April 2018, for sequences 


that were predicted to encode proteins similar to proinsulin. 


Initial searches used the tBLASTn algorithm (Gertz et al., 
2006) using previously characterized zebrafish (Danio rerio) 
proinsulin protein sequences (Irwin, 2004; Milewski et al., 1988) 
as queries. Subsequent tBLASTn, BLASTn, and BLASTp 
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searches used our putative proinsulin sequences as queries. 
Genome sequences of 55 fish (52 ray-finned (Superclass 
Actinopterygii), 2 cartilaginous (Superclass Chondrichthyes), 
and one lobe-finned (Class Sarcopterygii)) were examined, with 
52 of these contained in the NCBI Genome Data Viewer, and 
three (Atlantic cod, Gadus morhua; three-spined stickleback, 
Gasterosteus aculeatus; and the spotted green pufferfish, 
Tetraodon nigroviridis) only in the Ensembl database (9 were 
in both NCBI and Ensembl). A list of species with genomes 
used here, and their phylogenetic relationships, is presented 
in Figure 1. All sequences with E-scores below 0.001 were 
examined. As an initial step to assess homology and orthology, 
reciprocal BLASTp searches of the zebrafish proteome were 
conducted to examine the similarity of the putative proinsulin 
sequences. Sequences that were not more obviously similar to 
insulin-like growth factor-I (igff) or insulin-like growth factor-Il 
(igf2) were used for the following analyses. Here | suggest a 
standardized set of names for fish ins genes, based on the 
analysis described below, with insa, insb, and insc representing 
the three ins paralogs that originated early in teleost fish 
evolution, and insaa and insab representing duplicates of the 
insa gene shared by multiple species. Arabic numerals are 
used for recent lineage-specific duplicates and do not indicate 
necessarily indicate orthology. 


Phylogenetic analysis 

Prior to phylogenetic analysis, proinsulin (ins) coding 
sequences were aligned using the MAFFT algorithm 
(Katoh et al., 2001) as implemented at the Guidance2 
server (http://guidance.tau.ac.il/ver2/; Penn et al., 2010). 
Phylogenetic relationships were established using both 
maximum likelihood and Bayesian methods. PhyML 3.0 
(http://www.atgc-montpellier.fr/phyml/; Guindon et al., 2010) 
was used to generate maximum likelihood trees, with 500 
bootstrap replications, where the Smart Model Selection 
(SMS) (Lefort et al., 2017) option was used to identify the best 
fitting evolutionary model. Bayesian trees were constructed 
using MrBayes 3.2.6 (Huelsenbeck et al., 2001; Ronquist et 
al., 2012), with 5 000 000 generations and four simultaneous 
Metropolis-coupled Monte Carlo Markov chains sampled every 
100 generations under the same evolutionary model used for 
maximum likelihood. The first 25% of the trees were discarded as 
burn-in with the remaining samples used to generate consensus 
trees. Insulin coding sequences from cartilaginous (Callorhinchus 
milli (Elephant shark) and Rhincodon typus (Whale shark)) and 
lobe-finned fish (Latimeria chalumnae (Coelacanth)) were used 
as outgroups to root the trees for all insulin coding sequences. 
Subsequent analyses, using subsets of the sequence data, 
were conducted with similar approaches, and used sequences 
identified in our initial analyses as outgroups. 


Evolutionary analysis 

Genes adjacent to insulin genes in fish genomes were identified 
by examining genomic contigs in the NCBI and Ensembl 
databases using methods that | have previously used (Irwin, 
2002, 2012; Irwin & Mojsov, 2018). Briefly, the identity and 
orientation of genes neighboring insulin genes were identified 


from the genomic databases. 


The orthology of the genes 


adjacent to insulin genes was assessed (or confirmed if 
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Figure 1 Diversity and evolution of fish insulin genes 
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annotated) using BLASTp searches of the zebrafish and 
human proteomes (Kuraku & Meyer, 2012). 


Callorhinchus milii 
Rhincodon typus 
Latimeria chalumnae 
Lepisosteus oculatus 
Scleropages formosus 
Paramormyrops kingsleyae 
Clupea harengus 
Sinocyclocheillus anshuiensis 
Sinocyclocheillus grahami 
Sinocyclocheillus rhinocerous 
Cyprinus carpio 

Danio rerio 

Astyanax mexicanus 
Pygocentrus natterei 
Ictalurus punctatus 

Esox lucius 
Oncorhynchus kisutch 
Oncorhynchus mykiss 
Salmo salar 

Salvelinus alpinus 

Gadus morhua 
Hippocampus comes 
Boleophthalmus pectinirostris 
Monopterus albus 

Seriola dumerili 

Seriola lalandi 
Cynoglossus semilaevis 
Paralichthys olivaceus 
Lates calcarifer 
Acanthochromis polyacanthus 
Amphiprion ocellaris 
Stegastes partitus 
Oreochromis niloticus 
Haplochromis burtoni 
Neolamprologus brichardi 
Pygocentrus nyererei 
Maylandia zebra 

Poecilia formosa 

Poecilia. latipinna 
Poecilia mexicana 
Poecilia. reticulata 
Xiphophorus maculatus 
Fundulus heteroclitus 
Cyprinodon variegatus 
Nothobranchius furzeri 
Austrofundulus limnaeus 
Kryptolebias marmoratus 
Oryzias latipes 

Oryzias melastigma 
Labrus bergylta 
Larimichthys crocea 
Takifugu rubripes 
Tetraodon nigroviridis 
Noyothenia coriiceps 
Gasterosteus aculeatus 


I Cartilagenous fish 
Lobe-finned fish 


Ray-finned fish 


Teleosts 


Percomorpha 


Summary of the evolution of insulin genes in ray finned fish. Phylogenetic relationships of species with genome sequences is based on Betancur-R et al. (2017), 


with modifications in an area that is not confidently resolved to minimize the number of gene duplications or deletions. Some phylogenetic relationships are not 


resolved. Branches are not proportional to evolutionary time. Classification of fish is shown to the right. Solid diamonds indicate genome duplications. The 


numbers of different types of insulin gene found in each genome is listed, with numbers in brackets indicating incomplete sequences. Gene additions or gene 


losses are indicated by plus (+) and minus (—) before an a, b, or c, representing the insa, insb, and insc paralogs, respectively. The insa, insb, and insc paralogs 


originated near the genome duplication in the ancestor of teleosts and are indicated as +a, +b, and +c to refer to the two gene duplications. The duplication to 


generate the insaa and insab genes is indicated by +aa/ab. «ins indicated the gene duplication on the Lepisosteus oculatus lineage. «1pins, «ipa, and +1pb 


indicates retroposition events that generated processed pseudogenes. 


The strength of selective pressure acting upon the insulin 
coding sequences can be measured by comparing the relative 
rates of nonsynonymous (dy) and synonymous substitutions 
(ds). Sequences under stronger selective pressure for protein 
function will have lower nonsynonymous to synonymous 


(dw/ds) ratios. dw/ds rate ratios were obtained from analyses 
using RELAX (Wertheim et al., 2014) as implemented on the 


Datamonkey Adaptive Evolution server (http://datamonkey.org/; 
Weaver et al., 2018), which also tested for intensification 


or relaxation of the levels of selection on tested lineages. 
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Evidence for positive selection on branches of the phylogenetic 
tree was tested using aBSREL (Smith et al., 2015) on the 
Datamonkey Adaptive Evolution server (Weaver et al., 2018). 


Prediction of protein processing sites and alignment of 
peptide sequences 

Signal peptidase cleavage sites in the proinsulin protein 
sequences were predicted using SignalP 4.1 (http://www.cbs. 
dtu.dk/services/SignalP/; Petersen et al., 2011). | used two 
prediction programs to identify potential prohormone protease 
processing sites in the proinsulin sequences, NeuroPred (http: 
//neuroproteomics.scs.illinois.edu/neuropred.htm; Southey et 
al., 2006) and ProP 1.0 (http://www.cbs.dtu.dk/services/ProP/; 
Duckert et al., 2004), with "general PC prediction" selected. 
To visual the conservation of the processing sites and insulin 
peptide sequences, the proinsulin protein sequences were 
aligned using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/ 
clustalo/; Sievers et al., 2011). Consensus proinsulin A- 
and B-chain peptide sequences were generated from the 
alignments using WebLogo 3 (http://weblogo.threeplusone. 
com/; Crooks et al., 2004). 


RESULTS 


Numbers of insulin genes in fish genomes 

Searches of genome databases maintained by NCBI and 
Ensembl resulted in the identification of 168 insulin genes 
from the genomes of 55 species examined, with 141 of 
these predicted to be intact full-length coding sequences 
(Figure 1 and Supplementary Table S1). While only single copy 
insulin genes were found in the genomes of cartilaginous and 
lobe-finned fish, multiple genes were found in the genomes 
of most ray-finned fish. Within ray-finned fish, the number 
of identified insulin genes ranged from 1 (Sclerophages 
fromosus) to 8 (Sinocyclocheilus anshuiensis). As expected, 
fish species that experienced recent genome duplications, i.e., 
carp (Cyprinus carpio; Xu et al., 2014), members of the 
genus Sinocyclocheilus (Yang et al., 2016), and salmonids 
(Oncorhynchus kisutch, O. mykiss, Salmo salar, and Salvelinus 
alpinus; Lien et al., 2016) have the largest numbers of insulin 
genes. Phylogenetic analysis (see below) suggested the 
presence of three insulin paralogs in teleost fish: insa, insb, 
and insc. Of the 162 insulin genes identified in the 51 species 
of teleost fish (infraclass Teleostei), 92 were classified as insa 
(which includes insaa and insab genes), 57 as insb, and 13 
as insc, with 11, 12, and 3 of the sequences found to be 
incomplete, respectively, for these three paralogs (Figure 1 
and Supplementary Table S1). All 51 species of teleost fish 
possessed at least one /пѕа gene, with 45 having insb genes, 
and 10 with insc genes (Figure 1 and Supplementary Table 
S1). The 6 species of teleost fish where | failed to identify 
insb genes are distributed across the accepted phylogeny of 
fish (Figure 1; Betancur-R et al., 2017; Near et al., 2012), 
suggesting they were lost independently several times (or are 
missing from these genome assemblies). In contrast, almost 
all of the species lacking insc genes share a common ancestor 
(Figure 1, clade Euteleostei), suggesting a single gene loss event 
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explains the absence of insc genes in most teleost fish. The 
coding sequences of the insulin (ins) genes used in our analysis 
described below are presented in Supplementary Figure S1. 


Phylogeny of fish insulin genes 

To better understand the origin and evolution of the multiple 
insulin (ins) genes found in the genomes of ray-finned 
fish, their phylogenetic relationships, rooted using ins gene 
sequences from cartilaginous and lobe-finned fish, were 
inferred using maximum likelihood (Figure 2) and Bayesian 
(Supplementary Figure S2) methods. Both analyses generated 
similar phylogenies with independent duplications of the insulin 
gene on the spotted gar (Lepisosteus oculatus; infraorder 
Holostei) and teleost (infraorder Teleostei) fish lineages. While 
duplicated insa and insb genes were previously identified 
in teleost fish (Caruso & Sheridan, 2014; Hrytsenko et al., 
2016; Irwin, 2004), here | find that these species have a third 
insulin paralog, which | named insc. The species relationships 
within each paralog of insulin genes, given its limited power 
at resolving species relationships, is in general agreement 
with the accepted phylogeny of teleost fish (Betancur-R et al., 
2017; Near et al., 2012). This observation suggests that the 
triplication of the insulin gene occurred in an early teleost. The 
relationships among the insa, insb, and insc paralogs is poorly 
resolved in our phylogenetic analyses, with the maximum 
likelihood analysis suggesting a closer relationship between 
insa and insb, while the Bayesian analysis suggesting that insb 
and insc are most closely related. 

Most teleost species that possess insb and insc genes 
have a single copy of each gene, except those species 
that have experienced additional genome duplications (i.e., 
Cyprinus carpio, salmonids (O. mykiss, O. kisutch, Salmo salar, 
and Salvelinus aplinus), and the genera Sinocyclocheilus)). 
In contrast, a large number (30 of 51) of teleosts have 
multiple insa genes (Figure 1 and Supplementary Table S1). 
While some of these appear to be due to lineage-specific 
genome duplications (i.e., Cyprinus carpio, salmonids, and 
the genera Sinocyclocheilus) or single-gene duplications (e.g., 
Labrus bergylta and Larimichthys crocea), many are due to 
a duplication event that occurred on an early Percomorpha 
lineage (Figures 1, 2, and Supplementary Figure S2). The 
genomes of many Percomorpha fish contain two (or more) 
insa genes that form monophyletic clades that are largely 
in accord with the accepted species phylogeny (Betancur-R 
et al., 2017; Near et al., 2012). This observation suggests 
that a duplication of the insa gene, generating insaa and 
insab genes, occurred in an early lineage of Percomorpha 
fish (see Figure 1). In species that are inferred to be 
descendants of this Percomorpha-specific ins gene duplication, 
all but Nothobranchius furzeri possess an insaa gene, however, 
a larger number of species, such as tilapia (Oreochromis 
niloticus), some cichlids (e.g., Pundamilia nyererei and 
Maylandia zebra), and species of the genus Poecilia do not 
have an insab gene (Figures 1, 2 and Supplementary Figure 
S2 and Table S1). This suggests that the insab gene was lost in 
parallel in a number of species of Percomorpha (see Figure 1). 
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Figure 2 Phylogeny of fish insulin genes 
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Phylogenetic relationship of insulin coding sequences was inferred using maximum likelihood with the best fitting evolutional model GTR+G+I (GTR with gamma 


distribution and invariant sites) and bootstrapped with 500 replications. A similar phylogeny was generated using Bayesian methods (Supplementary Figure S2). 


Percent bootstrap support for lineages are shown. The tree in Newick format, with all branch lengths and support values, is presented in Supplementary Figure 


S3. The phylogeny was rooted with sequences from cartilaginous (Callorhinisus milli and Rhincodon typus) and lobe-finned (Latimeria chalumnae) fish. Full 


Species names are listed in Figure 1. insa, insb, and insc indicate the three clades of insulin genes found in teleost fish. 
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Duplications of fish insulin genes 


To aid in the assessment of the orthology of insulin genes, 
and to gain insight into the gene duplication events, | identified 
genes flanking the insulin genes in the genomic sequences of 
55 species of fish (Supplementary Table S2, examples given 
in Figure 3). A genome duplication occurred in the ancestor 
of teleost fish, the teleost fish-specific genome duplication 
(Glasauer & Neuhauss, 2014; Meyer & Van de Peer, 2005). 
Our phylogenetic analyses (Figure 2 and Supplementary S2) 
suggests that the three insulin paralogs found in teleost 
fish (insa, insb, and insc) originated on this lineage. I| 
duplicate genes originated through a genome duplication, 
then one might expect similar genomic neighborhoods for 
these genes (i.e., similar genes found adjacent to different 
ins genes), as observed for the duplicated proglucagon 
(gcg) and glucagon receptor (gcgr) genes of teleost fish 
(Irwin, 2014; Irwin & Mojsov, 2018). For the three insulin 
paralogs examined here, only very limited similarity was 
seen in the genomic neighborhoods (longer genomic regions 
were examined, with immediate neighborhoods summarized 
in Supplementary Table S2). A gene similar to tenm2 
(teneurin transmembrane protein 2) was found adjacent to insa 
genes in Exox lucius, O. kisutch, O. mykiss, Paramormyrops 
kingsleyae, Salmo salar, and Sclerophages formosus, and 
insb genes in Pygocentrus nattereri, Salvelinus alpinus, and 
Sinocyclocheilus rhinocerous (Supplementary Table S2). This 
might suggest that the insa and insb genes originated via 
the teleost-fish specific genome duplication. However, as the 
linkage between the tenm2 and ins genes was only observed 
in a small number of species, this would require large number of 
independent losses (or genomic rearrangements) of tenm2-like 
genes to yield the absence of a linked tenm2-like with most 
ins genes, therefore the evidence for the origin of insa and 
insb by a genome duplication event is not strong. No genes 
near insc were similar to those adjacent to insa or insb, thus 
it more likely originated as a single gene duplication event 
rather than through a genome duplication. If the insa and 
insb genes originated via a genome duplication, an addition 
single-gene duplication would still be needed to explain the origin 
of insc. The gene neighborhood data does not help resolve 
the relationships among the three ins genes. Intriguingly, the 
genomic neighborhoods of none of the insulin gene paralogs 
in teleost fish is similar to the typical vertebrate insulin gene 
neighborhood, where insulin is flanked by insulin-like growth 
factor 2 (igf2) and tyrosine hydroxylase (th) genes (Patton et al., 
1998), although similar neighborhoods were found for the insulin 
genes in cartilaginous and lobe-finned fish (Supplementary Table 
S2). This suggests that the insulin genomic neighborhood was 
rearranged on the ancestral ray-finned fish lineage. 


Orthologous genes often have similar genomic 
neighborhoods, thus can these neighborhoods can be used 
to assess orthology when gene sequence similarity is limited 
(Kurokawa et al., 2005). When the genomic neighborhoods 
surrounding insb and insc genes were examined (Figure 3 
and Supplementary Table S2) shared features were identified 
for each gene. All insb genes are adjacent to a hmmr 
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(hyaluronan mediated motility receptor)-like gene, while all 
insc genes are flanked by prr12 (proline rich 12) and scaf1 
(SR-related CTD associated factor 1)-like genes (as shown 
for the Clupea harengus insb and Danio rerio insc genes 
in Figure 3), suggesting that their genomic neighborhoods 
are largely conserved since they originated in the early 
teleost. In contrast, several different genomic neighborhoods 
were observed for the insa genes (Supplementary Table S2, 
examples shown in Figure 3). The most phylogenetically 
widespread genomic neighborhood has a nipal (magnesium 
transporter NIPA3-like)-like gene adjacent to the insa gene 
(Supplementary Table S2) and is found in species representing 
Osteoglossomorpha (e.g., Scelerophages formosus, shown 
in Figure 3), Protacanthopterygii (e.g. Esox lucius and 
Salmo salar), Paracanthopterygii (e.g., Gadus morhua), and 
Percomorpha (e.g., Takifugu rubripes and Oryzias latipes). 
The only major group of teleost fish that does not have an insa 
gene adjacent to nipal are those of Otomorpha (e.g., Danio 
rerio, shown in Figure 3), suggesting a change in the gene 
order occurred in this group and that the linkage of insa with 
nipal is ancestral. Within Percomorpha (Supplementary Table 
S2), the insaa genes are found adjacent to nipal-like genes (as 
in Scelerophages formosus, Figure 3), while the insab genes 
are near a sh3pxd2b (SH3 and PX domains 2B)-like genes 
(e.g., Haplochromis burtoni, Figure 3). The linkage of the insab 
genes with sh3pxd2b strengthens the conclusion that all insab 
genes are orthologous and indicates that insaa is located at the 
ancestral genomic location and was the source (locus-of-origin) 
for the insab gene that was inserted into a new genomic locus. 
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Figure 3 Genomic neighborhoods surrounding fish insulin 
genes 

Examples of genomic neighborhoods found for insulin paralogs in ray-finned 
fish. Genomic neighborhoods are listed in Supplementary Table S2. A: 
Typical genomic neighborhood for insa and insaa genes as illustrated from 
Sceleropages formosus. B: Genomic neighborhood for insa genes in Danio 
rerio and close relatives. C: Typical genomic neighborhood for insab genes 
as illustrated from Haplochromis burtoni. D: Typical genomic neighborhood 
for insb genes as illustrated from Clupea harengus. E: Typical genomic 
neighborhood for insc genes as illustrated from D. rerio. See Supplementary 
Table S2 for full gene names. 


Lineage-specific changes in insulin gene number 


The gene duplication events described above explain most, 
but not all, of the insulin genes identified (Figure 1 and 
Supplementary Table S1). Many of the additional genes 
can be explained by tandem gene duplication, which lead to 
two insulin genes arranged head-to-tail in the genome (see 
Supplementary Tables S1, S2) and includes the insat and 
insa2 genes in Parmormyrops kingsleyae and in Larmichthys 
crocea, the іпѕааї and insaa2 genes in Hippocampus 
comes, and the insab1 and insab2 genes in Boleopththalmus 
pectinirostris. Both Cyprinus carpio and Sinocyclocheilus 
rhinocerous have tandemly arranged insb genes (insb3 and 
insb4 in C. carpio and insb2 and insb3 in S. rhinocerous, 
Supplementary Table S2) that might have originated prior to 
the divergence of these two species (although, they would 
then have to be lost in the two other Sinocyclocheilus 
species). Sinocyclocheilus anshuiensis has two pairs of 
tandemly arranged genes, insa1/insa2 and insa3/insa4, which 
either originated prior to this species genome duplication (and 
then were lost in the closely related species that share this 
genome duplication), had parallel tandem duplications, or is an 
assembly error (the tandemly arranged genes are identical in 
sequence). The ins? and ins2 genes of Lepisosteus oculatus 
might also have originated via a tandem duplication, which was 
then followed by rearrangements so that these lineage-specific 
gene duplicates are now found in opposite orientations on 


either side of a hmmr-like gene (see Supplementary Table S2). 


More complex duplication events, resulting in genes moving to 
new locations are need to explain the remaining genes, with 
some being extremely recent (e.g., the identical insaaí and 
insaa2 genes in Labrus bergylta) or more ancient (e.g., insat 
and insa2 gene of Esox lucius). 

While duplication of genomic DNA is a frequent mechanism 
for the origin of duplicate genes, they can also be generated 
by retroposition of cDNA generated from mRNA transcripts 
(Long et al., 2013). Indeed, the well-characterized duplicate 
insí gene in mice was generated from an incompletely 
processed mRNA transcript (Shiao et al., 2008; Soares et 
al., 1985). Typically, insulin genes contain two coding exons 
that must be spliced together to create an intact coding 
sequence. Here | identified 6 insulin genes that are likely 
generated by reverse-transcribed mRNAs, as these genomic 
sequences do not have an intron interrupting the coding 


sequence (see Supplementary Figure S4 for an example). 


The 6 retro-processed insulin genes can be explained by 
three integration events: (1) between ankrd6 (ankyrin repeat 
domain 6) and /yrm2 (LYR motif containing 2)-like genes on the 
Lepisosteus oculatus lineage to generate wins; (2) within an 
apbb3 (amyloid beta precursor protein binding family 8 member 
3)-like gene on the Monopterus albus lineage to generate 
wWinsb; and (3) within a psd3 (PH and SEC7 domain containing 
protein 3)-like gene on the ancestral lineage for salmonids 
(Oncorhynchus mykiss, O. kisutch, Salmo salar, and Salvelinus 
aplinus) and the Northern pike (Esox lucius) that generated 
winsa іп O. kisutch, S. salar, and S. alpinus (but was not 
found in O. mykiss) and an intact coding sequence (insa2) in 


Esox lucius (Figure 1 and Supplementary Figure S4 and Table 
S2). The maintenance of an intact open reading frame in Exox 
lucius insa2 (Supplementary Figure S4) suggests that this 
retropositioned sequence, like the mouse ins1 gene (Soares et 
al., 1985), is still functional. However, this sequence has lost its 
coding potential in salmonid fish and is a pseudogene (Figure 1 
and Supplementary Table S1). Additional processed genes 
might also exist, as several incomplete genes (e.g., insb4 from 
Cyprinus carpio and insa3 from Paramormyrops kingsleyae) 
were identified (see Supplementary Table S1) that were located 
at unique location of the genome (see Supplementary Table 
S2) consistent with being inserted processed cDNAs, but the 
sequences were similar to either exon 1 or 2 (and not both), 
thus cannot be distinguished from a sequence generated by an 
DNA-mediated incomplete gene duplication event. 


Gene for insulin a (insa) is evolving under greater 
evolutionary constraint 


Our genomic and phylogenetic analyses of fish genome 
sequences demonstrated that teleost fish have three 
paralogous insulin genes. In contrast, most other vertebrate 
species only have one (Chan & Steiner, 2000; Conlon, 
2000b). The increased number of insulin genes in teleost 
fish raises the possibility that they have been diversified 
to: (1) subfunctionalize distinct functions of insulin, (2) 
neofunctionalize to acquire novel functions, or (3) lose function 
(pseudogenize). As a first step to explore possible changes in 
the biological roles of these distinct insulin genes | assessed 
the selective pressures acting upon the sequences. This 
can be done by comparing the rates of nonsynonymous 
(dy) to synonymous (ds) substitutions (Yang & Bielawski, 
2000). If the protein encoded by a gene has lost its function, 
then it would be expected that there would be no selection 
against nonsynonymous substitutions, and that rates of 
nonsynonymous and synonymous substitutions would be the 
same rate (і.е., dw/dg-1, neutral evolution). If the protein 
encoded by these genes still had a function, then selection 
should act against a subset of nonsynonymous substitutions 
(the deleterious mutations) and be lower than the synonymous 
rate (i.e., dy/ds«1, purifying selection). Rarely, one might see 
dw/ds»1, which would indicate positive selection for change 
in amino acid sequence (Yang & Bielawski, 2000). If protein 
coding sequences of genes are being maintained for different 
functions, with different parts of the sequence being important 
for these functions, then one might see differing levels of 
selective constraint. 


| used the program RELAX (Wertheim et al., 2014) to 
assess the levels of selective pressure (dw/ds) for each of 
the three insulin gene paralogs, and to determine whether 
the selective pressure was intensified or relaxed compared to 
the other insulin genes, with sequences from a non-teleost 
ray-finned fish, Lepisosteus oculatus, used as outgroup. When 
all sequences were analyzed, insa genes showed the stronger 
levels of purifying selection (dy/ds=0.2666, 0.3332, and 
0.3258 for insa, insb, and insc, respectively; Table 1). If 
| restricted these analyses to either species that have all 
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three paralogs, or species that only have a single copy 
of these paralogs, to attempt to minimize species-specific 
effects (changes seen in species that only have some of the 
genes) or due to differences in gene number, the difference 
in the selective constraints acting on insa vs. insb and insc 
became even more pronounced (Table 1). The results from 
RELAX also indicated that selection intensification, compared 
to the other insulin sequences, occurred on insa, while 
selection relaxation occurred for insb (Table 1). No significant 
change in selection intensity was seen for insc vs. all other 
insulin sequences. Similar results, with insa showing greater 
constraint, were seen in comparisons of insa and insb genes 
in species that had both of these genes (Table 1). These 
results suggest that the insa paralog is under greater purifying 
selective constraint than either the insb or insc paralog, but 
also importantly demonstrate that both insb and insc are 
continuing to experience purifying selection. Thus, it can 
be concluded that the protein coding sequence of insa is 
under the greatest level of selective constraint, however, the 
coding sequences of both insb and insc are also experiencing 
selective constraints consistent with their protein products 
having essential biological functions. 


Table 1 Differences in the selective pressure (dy/ds) acting on 
insulin paralogs in teleosts 


Sequences tested? insa insb insc 





All 0.26665 (80) 0.3332° (44) 0.3258 (10) 
Species with insa, 


: : 0.2120° (11) 0.3309° (9) 0.3406 (9) 
insb, and insc 
А ith sinal 
Species with single copy O 1944 (4) 0.3164 (4) — 03129 (4) 
insa, insb, and insc 
Species with single — 1664 (17) 0.3171° (17) N/A 
copy insa and insb 
Species with single 
copy insa and 0.1366? (13)  0.3191* (13) N/A 


insb, but no insc 


Numbers in brackets is the numbers of sequences used. 2: The 
two complete Lepisosteus acaulatus ins sequences were used as 
outgroups. b: Test for selection intensification was significant. °: 
Test for selection relaxation was significant. N/A: Not applicable. 


Adaptive evolution of insulin genes 

While the protein products of the three insulin genes are 
being maintained for function, they might not be the same 
function. To examine whether any of the sequences might 
have gained new functions | tested for evidence of positive 


selection on lineages using aBSREL (Smith et al., 2015). 


Only one lineage showed significant evidence for episodic 


diversifying selection, the ancestral lineage for the insc genes. 


Inspection of the phylogenetic tree generated from all of the 
insulin coding sequences suggest that the ancestral branch 
for the insc genes is longer than for the insa and insb genes 
(Figure 2 and Supplementary Figure S2), consistent with a 
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greater number of nonsynonymous substitutions driven by 
positive selection. No evidence for positive selection was found 
on the ancestral lineages for the insa or insb genes. These 
results might suggest that insc has acquired a new function 
(neofunctionalization) driven by positive selection, while insa 
and insb have been retained due to subfunctionalization of 
insulin functions between the two genes. 


Relaxed selection on the insab duplicate of the insulin a 
(insa) gene 

An intriguing observation from our analysis of the selective 
constrains acting upon insa, insb, and insc paralogs was that 
the calculated value for the selective constraint acting on insa 
varied considerably depending upon which sequences were 
used for the analysis (Table 1). The ratio of the nonsynonymous 
to synonymous rates was higher when all insa sequences were 
used than if sequences were used only from species that had 
all three paralogs or had only single copies of the insa and 
insb paralogs (Table 1). These observations suggest that 
inclusion of insa sequences from species that have multiple 
insa sequences yields higher estimates of the dy/ds ratio due 
to some of these sequences having lower levels of sequence 
constraint. As the duplication to generate the insaa and insab 
genes is a major source of the multiple insa genes in teleost 
fish (see Figures 1, 2 and Supplementary Figure S2 and Table 
51), | compared the selective constraints acting upon insaa 
and insab genes (Table 2). The insaa coding sequences 
were found to be under selective constrains (Table 2) similar 
to those of other insa sequences, and especially those 
from species that had single copy insa, while the insab 
sequences displayed the lowest levels of constrains seen for 
any fish ins sequence (Tables 1, 2). Test for intensification 
or relaxation of selective constraint showed that insaa was 
under significant intensification of selective constraint, while 
insab was significant relaxation of selective constraint was 
demonstrated for insab (Table 2). These results suggest that 
insaa retains the function of insa, which would be consistent 
with it being at the locus of origin, while insab is evolving with far 
less constraint, to the extent that it has been lost on a number of 
lineages (e.g., Oreochromis niloticus, Pundamilia nyererei, M. 
zebra, and species of the genus Poecilia, see Figure 1). This 
might suggest that insab is not essential. 


Table 2 Differences in the selective pressure (dw/ds) acting on 
insaa and insab genes 


Sequences tested? insaa insbb 





All insa 0.1691 (18) 0.4766° (14) 


Species with one 


0.12510 (9) 
copy of each 


0.5139° (9) 





Numbers in brackets is the numbers of sequences used. ?: The 
two complete Lepisosteus acaulatus ins sequences were used as 
outgroups. ^: Test for selection intensification was significant. °: 
Test for selection relaxation was significant. 


Processing of proinsulin sequences 

To be functional, the protein product encoded by insulin 
genes need to be proteolytically processed and secreted to 
generate two-peptide chain insulin molecules (Steiner et al., 
1996, 2009). Our analysis was initially only focused on 
sequences that showed similarity to previously characterized 
proinsulin sequences and had intact coding sequences, i.e., 
had an initiation codon, a termination codon and intact open 
reading frame with greater similarity to proinsulin than to other 
insulin-like sequences (e.g., insulin-like growth factors). To 
determine whether the encoded protein sequences could be 
secreted and properly processed | searched for potential signal 
peptidases and proprotein processing sites using programs 
that predict these sites (Duckert et al., 2004; Petersen et al., 
2011; Southey et al., 2006). 

The coding sequences of intact proinsulin open reading 
frames from all ray-finned fish identified here are predicted to 
have functional signal peptides (Table 3 and Supplementary 
Table S3 and Figure S6), thus should be able to be secreted. In 
addition, all of the proinsulin protein sequences encoded by the 
insa, insab, insb, and insc genes have potential prohormone 
protease cleavage sites that could yield typical two-chain 
insulin hormone molecules (Table 3 and Supplementary Table 
S3 and Figure S6). In contrast, however, only one of the 14 
proteins encoded by insab genes is predicted to potentially 
produce a two-chain insulin (Table 3). For the insab sequences, 
only 2 have potential B-chain/C-peptide processing site and 9 
have potential C-peptide/A-chain processing sites (Table 3 and 
Supplementary Table S3 and Figure S6), indicating that neither 
site is conserved. Both NeuroPred (Southey et al., 2006) and 
ProP (Duckert et al., 2004) predicted similar sites, but ProP 
also provided a list of other potential processing sites that did 
not score well enough to be predicted sites (see Supplementary 
Table S3). The alternative potential sites identified by ProP in the 
insab proteins still would not generate 2-chain insulin molecules 
similar to functionally characterized insulin molecules (Steiner et 
al., 1996, 2009) as they would generate very short B-chains, 
often missing residues essential for function. 


Table 3 Secretion and processing of fish proinsulin proteins 
Signal B-C C-A 





2-chain 





Geng peptide? junction? junction? insulin° 
insa 35/358 35/35 35/35 35/35 
іпѕаа 32/32 32/32 32/32 32/32 
insab 14/14 2/14 9/14 1/14 
insb 45/45 44/45 45/45 44/45 
insc 10/10 10/10 10/10 10/10 





a: Sequences with signal peptidase processing sites. ^: Sequences 
with proprotein processing sites for the B-chain/C-peptide and 
C-peptide/A-chain processing sites. ^: Number that can be secreted 
and properly processed to yield two-chain polypeptides. 3: Number 
with the site/total number of sequences. 


The hormone insulin is a two-chain peptide molecule held 


together by disulfide bridges (Steiner et al., 1996, 2009). 
To confirm that two-chain molecules could be produced 
from the proinsulin protein sequences | generated consensus 
sequences for the A- and B-chains for the different types 
of insulins (Figure 4). The lengths of the predicted A- and 
B-chains were similar between the different types of insulins, 
with more variation, at both the N- and C-termini, within proteins 
encoded by a type of gene than between types of genes 
(Figure 4 and Supplementary Figure S6). When cysteine 
residues were examined one of the two insc sequences from 
Sinocyclocheilus grahami (insc2) had a cystine replaced by an 
arginine residue. It is possible that the S. graham insc2 gene 
does not encode a functional insulin, but this is compensated 
by the presence of an insc1 gene that encodes an insulin that 
can be secreted processed and has all six cysteine residues. 
Surprisingly, when | examined A- and B-chain sequences 
homologous to those of other insulins in the protein encoded 
by insab genes, all of the cysteine residues were perfectly 
conserved, however much of the rest of the sequence showed 
greater variation than within other types of insulins (Figure 4 
and Supplementary Figure S6). This raises the possibility that 
insab could still fold in a similar fashion as typical insulin, thus 
be transported to insulin secretory granules, but then not be 
processed into an active hormone (Liu et al., 2015). 


DISCUSSION 


Searches of diverse fish genomes, including in the extensively 
characterized zebrafish (Danio rerio) genome, unexpectedly 
revealed the existence of three insulin gene paralogs (Figure 1 
and Supplementary Figure S1 and Table S1). Previously, two 
paralogous insulin genes, insa and insb, had been found in 
the zebrafish genome that have distinct expression patterns 
(Irwin, 2004; Papasani et al., 2006). Here | found insa genes 
in the genomes of all teleost fish (Figure 1 and Supplementary 
Figure S1 and Table S1), which encode the previously isolated 
insulin hormone sequences (Caruso & Sheridan, 2011; Conlon, 
2000b) that regulate glucose metabolism. Given the central 
role of insulin in the regulation of metabolism (Caruso & 
Sheridan, 2011; Róder et al., 2016; Weiss, 2009), it is not 
surprising that this gene was found in all fish. The presence of a 
high selective constraint (Tables 1, 2) on the protein sequences 
encoded by insa and insaa genes supports the conclusion 
that this sequence is essential. A surprising discovery was 
the presence of a duplicate of the insa gene, the insaa and 
insab genes, in a large number of teleost fish (Figure 1 
and Supplementary Table S1). The protein coding sequence 
encoded by insab genes have the highest nonsynonymous 
substitution rate (Table 2), and thus least selective constraint, 
yet the consensus sequences of the potential A- and B-chains 
of insulin predicted from these genes have perfectly conserved 
cystine residues (Figure 4). This suggests that these genes 
might encode proinsulin-like proteins that could properly fold 
but would not be properly processed in secretory granules to 
produce active insulin molecules (Liu et al., 2015). Thus, insab 
proteins might compete with insaa proteins for processing 
enzymes and modulate the release of functional insulin 
molecules in these species. 
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Figure 4 Consensus fish insulin peptide sequences 

Consensus sequences for the A-chain (A) and B-chain peptides (B) for insulin encoded by insa, insaa, insab, insb, and insc genes. Consensus sequences were 
generated from aligned predicted protein sequences (see Supplementary Figure S5) of genes with complete open reading frames (Larimichthys crocea insa2 
was removed as it may not be functional). Hight of the residue is proportional to representation among the compared sequences, while residue width represents 


proportion of sequences without a gap (thin residues are gaps in many sequences). 


The insb gene was found in the genomes of almost all teleost Table S1). The selective forces acting on insb are weaker 
fish species examined, demonstrating that it is conserved than those acting on insa (Table 1) suggesting that these 
in these fish (Figure 1 and Supplementary Figure S1 and two genes have different functions. Previous work indicates 
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that insb is not appreciably expressed in the pancreas, and 
instead is more predominantly expressed in early development 


(Hrytsenko et al., 2016; Irwin, 2004; Papasani et al., 2006). 


The analysis of selective constraints and these expression 
results are consistent with speculation that insb is not primarily 
involved in the regulation of blood glucose level, but instead has 
a role in development (Hrytsenko et al., 2016; Papasani et al., 
2006). Given the conservation of the insb gene across teleost 
fish, it likely acquired this role soon after its origin and has been 
conserved. 

Unexpectedly | found an insc paralog, and it was present in 
the well characterized zebrafish (D. rerio) genome (Figure 1 
and Supplementary Figure S1 and Table S1). Upon closer 
examination, the insc gene in both the Ensembl (gene ID: 
ENSDARGO00000096862) and NCBI (GenBank accession No.: 
100534937) zebrafish genomes fail to predict an intact open 
reading frame, with the intact insc coding sequence only 
coming from the linked mRNA sequence (XM 009299616.3) in 
the NCBI database (see Supplementary Table S1). | could find 
no raw sequence data (e.g., ESTs) that supports the existence 
of an intact open reading frame for zebrafish insc, as all mRNA 
derived sequences (GenBank accession Nos. EH533807.1 
and EH507658.1, forward and reverse sequences of the same 
EST clone) also contained the deletion found in the genomic 
sequence. This suggests that the zebrafish insc gene does 
not encode an intact protein product. However, our analysis 
of the selective forces acting upon the insc coding sequences 
indicate that the nonsynonymous rate of substitutions (dy) is 
much lower than the rate of synonymous substitutions (ds) 
consistent with selection to maintain a protein coding sequence 


(Table 1) in most, if not all, species that have intact insc genes. 


The insc gene, in contrast to insa and insb, is restricted to the 
genomes of only a few lineages, with this gene been lost from 
the genome of most teleost fish (Figure 1 and Supplementary 
Table S1). This distribution of the gene among fish suggests 
that it does not encode a function that is essential in most 
Species, and its loss might have been compensated by the 
presence of paralogs of insulin in the fish genomes. Expression 
of insc, based on a single EST clone (with sequences from both 
ends) that was identified in zebrafish, is found in the “gut and 
internal organs" of adults, thus potentially it could be expressed 
in the pancreas and be replaced by insa. Further studies on the 
expression of insc in other fish are needed to identify its role in 


physiology. 
CONCLUSION 


Searches of ray-finned fish genomes have demonstrated that 
they contain more insulin genes than previously appreciated 


and suggest that the roles of these genes have diversified. 


The diversification of the functions of insulin parallels the 
diversification of the proglucagon-encoded peptides found in 
fish (Irwin & Mojsov, 2018). Similarly, parallel changes in 
the biological functions of insulin and proglucagon-derived 
peptides have been observed in hystricomorph rodents (i.e., 
guinea pig and relatives), where changes in sequence and 
action of insulin have compensating changes in glucagon 


(Seino et al., 1988). 
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ABSTRACT 


Rana kunyuensis is a species of brown frog that lives 
exclusively on Kunyu Mountain, Yantai, China. In the 
current study, a 279-bp cDNA sequence encoding a novel 
antimicrobial peptide (AMP), designated as amurin-9KY, 
was cloned from synthesized double-strand skin cDNA 
of R. kunyuensis. The amurin-9KY precursor was 
composed of 62 amino acid (aa) residues, whereas the 
mature peptide was composed of 14 aa and contained 
two cysteines forming a C-terminal heptapeptide ring 
(Rana box domain) and an amidated C-terminus. These 
structural characters represent a novel amphibian AMP 
family. Although amurin-9KY exhibited high similarity to 
the already identified amurin-9AM from R. amurensis, 
little is known about the structures and activities of 
amurin-9 family AMPs so far. Therefore, amurin-9KY and 
its three derivatives (amurin-9KY1—3) were designed and 
synthesized. The structures and activities were examined 
to evaluate the influence of C-terminal amidation and 
the heptapeptide ring on the activities and structure of 
amurin-9KY. Results indicated that C-terminal amidation 
was essential for antimicrobial activity, whereas both 
C-terminal amidation and the heptapeptide ring played 
roles in the low hemolytic activity. Circular dichroism 
(CD) spectra showed that the four peptides adopted an 
a-helical conformation in THF/H2O (v/v 1:1) solution, 
but a random coil in aqueous solution. Elimination 
of the C-terminal heptapeptide ring generated two 
free cysteine residues with unpaired thiol groups, 


198 Science Press 


which greatly increased the concentration-dependent 
anti-oxidant activity. Scanning electron microscopy (SEM) 
was also performed to determine the possible bactericidal 
mechanisms. 


Keywords: Antimicrobial peptides; Rana kunyuensis; 
Amurin-9KY; Heptapeptide ring; C-terminal amidation; 
Structure activity relationship 


INTRODUCTION 


Antimicrobial peptides (AMPs), which are small, cationic and 
amphipathic peptides widely distributed throughout organisms, 
are evolutionarily ancient weapons against environmental 
pathogens (Zasloff, 2002). They are important components of 
innate immune systems and play key roles in the anti-infective 
immune responses of organisms (Radek & Gallo, 2007). 
According to previous studies, AMPs possess strong and diverse 
antimicrobial activities against bacteria, fungi, viruses and even 
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protozoa. The activities and specificities of AMPs can be affected 
by various factors such as molecular weight, sequence, charge, 
conformation, hydrophobicity and amphipathicity (Brogden, 
2005). The antimicrobial mechanisms are diverse among 
different AMPs. However, most AMPs are thought to function 
by forming pores in the membranes of target microorganisms, 
ultimately leading to disruption of cellular integrity (Nicolas, 2009). 

Frogs belonging to the family Ranidae are widely distributed 
around the world, except for the Polar Regions, southern 
South America and most of Australia (Conlon et al., 2004). 
So far, 347 ranid species have been identified worldwide, 
with most living in moist environments surrounded by diverse 
pathogens. As a result, they are continuously threatened 
by pathogenic invasions and have therefore evolved effective 
immune systems (Conlon et al., 2004; Zasloff, 2002). Among 
these systems, AMPs play key roles in the anti-infective 
immune response of ranid frogs (Kreil, 1994). То date, 
hundreds of AMPs have been characterized from various ranid 
species. According to their primary structures, they are divided 
into dozens of different families, including gaegurins, brevinins 
(1 and 2), ranalexins, ranatuerins (1 and 2), esculentins (1 
and 2), palustrins, japonicins (1 and 2), nigrocin-2, tigerinins, 
temporins and odorranains (A-W) (Conlon, 2004; Duda et al., 
2002; He et al., 2012; Li et al., 2007; Matutte et al., 2000). 

Brown frogs, also known as wood frogs, are a group of 
ranid species belonging to Hana. Currently, there are 14 
recognized brown frogs in China distributed in 30 provinces. 
Rana kunyuensis, an endemic brown frog species in China, 
exclusively lives on Kunyu Mountain, Yantai, Shandong 
Province. In the present study, we identified a novel AMP family 
member, named amurin-9KY, from the skin secretions of R. 
kunyuensis. The structure-function relationship of amurin-9KY 
was clarified to establish the effect of C-terminal amidation and 
the heptapeptide ring on the biological activities and advanced 
structures of amphibian AMPs. 


MATERIALS AND METHODS 


Frog collection 

Adult specimens of R. kunyuensis (n=5, weight range 5—10 g) 
were captured on Kunyu Mountain, Yantai, Shandong Province, 
China. The frogs were housed in plastic box and fed with yellow 
mealworms in the lab for one week until the experiment. 


Total RNA extraction, cDNA synthesis and amurin-9KY- 
encoding cDNA screening 

An individual frog (female, 8 g) was washed with distilled water 
twice and anaesthetized with ice before being euthanized with a 
needle. The dorsal skin was quickly removed with a small pair 
of scissors and placed into a pre-cooling homogenizer. The 
skin was then immediately homogenized into powder in liquid 
nitrogen and mixed with Trizol reagent (Invitrogen, CA, USA). 
The subsequent procedures were carried out according to the 
manufacturer's instructions and the extracted total RNA was 
preserved in liquid nitrogen until use. All animal experimental 
protocols in the present study were approved by the Animal 
Care and Use Ethics Committee of Soochow University. 


Double-strand cDNA was synthesized using an In-Fusion 
SMARTer™ Directional cDNA Library Construction Kit 
(Clontech, USA). The primers used for first strand synthesis 
were 3' In-Fusion SMARTer CDS Primer, 5'-CGGGGTACG 
ATGAGACACCA(T);9 VN-3’ (N=A, C, G, or T; V=A, G, or C) 
and SMARTer V Oligonucleotide, 5’-AAGCAGTGGTATCAA 
CGCAGAGTACXXXXX-3’ (X=undisclosed base іп the 
proprietary SMARTer oligo sequence). The enzyme used 
for first strand synthesis was SMARTScribe™ Reverse 
Transcriptase and was supplied by the kit. The second 
strand was amplified by 50x Advantage 2 Polymerase Mix 
using 5’ PCR Primer Il A and 3’ In-Fusion SMARTer PCR 
Primer, 5’-CGGGGTACGATGAGACACCA-3’. The synthesized 
double-strand cDNA was stored at —80 °C until use. 

The synthesized double-strand cDNA was used as a 
template to screen the cDNAs encoding the AMPs of R. 
kunyuensis. A sense degenerate oligonucleotide primer 
(RanaAMP, 5’-CCAAAGATGTTSMCCWYGAAG-3’, M=A or C; 
W-=A or T; Y=C or T), designed according to the conserved 
signal peptide domain sequences of previously characterized 
AMPs from the skin of ranid frogs, and a 3’-antisense primer 
(3’-PCR, 5’-CGGGGTACGATGAGACACCAT-3’) were used for 
PCR analysis. The PCR procedure was as follows: 4 min of 
denaturation at 95 °C; 30 cycles: denaturation at 95 °C for 
30 s, annealing at 57 °C for 30 s, extension at 72 °C for 1 
min; final extension at 72 °C for 10 min. The PCR product was 
purified by agarose gel electrophoresis, cloned into pMD™19-T 
vector (Takara, Dalian, China) and transformed into E. coli for 
sequencing. 


Peptide synthesis 

Amurin-9KY and its derivatives (amurin-9KY1-3) were 
chemically synthesized on an Applied Biosystems model 433A 
peptide synthesizer (ABI, USA) according to the manufacturer's 
standard protocols. The crude peptides were purified by 
reversed-phase high performance liquid chromatography 
(RP-HPLC) to a final purity higher than 95% and their identities 
were confirmed by automated Edman degradation and 
matrix-assisted laser desorption/ionization—time-of-flight mass 
spectrometry (MALDI-TOF-MS). The intra-peptide disulfide 
bridge and N-terminal amidation were accurately formed in the 
synthesis process. 


Antimicrobial assay 

In total, eight microbial strains, including gram-positive 
bacteria, gram-negative bacteria and fungi, were used for 
the antimicrobial assay. The assays included inhibition 
zone examination and minimum inhibitory concentration (MIC) 
determination and were conducted as described in our previous 
papers (He et al., 2012; Lu et al., 2010). Briefly, microbes were 
inoculated in Mueller-Hinton broth (MH) and incubated at 37 °C 
to the exponential phase. The inoculum was then diluted with 
fresh MH broth to 10€ CFU/mL, and 50-uL bacterial dilutions 
were mixed with serial dilutions of peptides (50 uL) in 96-well 
microtiter plates. The plates were incubated at 37 °C for 18 h 
and the minimum concentrations at which no visible growth of 
bacteria occurred were recorded as MIC values. 
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Hemolytic assay 

Hemolytic assay was conducted as previously reported (Wang 
et al., 2011). Fresh human erythrocytes were collected from a 
healthy donor and washed with 0.9% saline. The erythrocytes 
were re-suspended with 0.9% saline and incubated with serial 
peptide dilutions at 37 °C for 30 min. The mixtures were then 
centrifuged at 2 000 r/min at room temperature for 5 min, after 
which the supernatants were removed and their absorbance 
was measured at 540 nm. We used 1% (v/v) Triton X-100 to 
determine the 100% hemolysis, with 0.9% saline used as the 
negative control. 


Anti-oxidant assay 
As a stable nitrogen radical, 2, 2-diphenyl-1-picrylhydrazyl 
(DPPH, Sigma, USA) has been widely used to examine the 


anti-oxidant activities of biological samples, drugs and foods. 


In the present study, the DPPH radical scavenging activities 
of the peptides were determined according to our previously 
described method (Zhang et al., 2012). Briefly, 92 рі of 
6x10-5 mol/L DPPH radical dissolved in methanol was mixed 
with 8 uL of peptide solution (2 mg/mL, mass ratio of peptide 
to DPPH of ~3:1). The mixture was kept in the dark for 30 min 
at room temperature and then the amount of reduced DPPH 
was quantified by measuring the decrease in absorbance at 
520 nm. The DPPH radical scavenging percentage (5%) 
of the peptides was calculated according to the formula: 
S%=(Abiank—Asample) x 100/Apiank- Deionized water was used as 
the negative control. 


Circular dichroism (CD) spectroscopy 
To investigate the secondary structure of amurin-9KY and its 
derivatives (amurin-9KY1—3), CD spectroscopy was performed 


using a Jasco J-810 spectrophotometer (JASCO, Japan). 


Samples with a constant peptide concentration of 0.12 mg/mL 
were prepared in two different solvents, i.e., water and 5096 
(v/v) trifluoroethanol (TFE)-water, and added to a quartz optical 
cell with a path length of 0.1 cm at 25 °C. The spectra were 
averaged over three consecutive scans, followed by subtraction 
of the CD signal of the solvent. 


Scanning electron microscopy (SEM) 

The surface morphologies of the AMP-treated microbes were 
observed by SEM, which can partly reveal the antimicrobial 
mechanisms of AMPs. In the present study, S. aureus 
ATCC25923 was used to evaluate the potential antimicrobial 
mechanism of amurin-9KY. The experiment was carried out 


according to the method described by Lu et al. (2010). 


The concentration of amurin-9KY used in the experiment 
was 1xMIC. The treatment conditions of the bacterial and 
amurin-9KY mixture were 37 °C for 30 min. After that, samples 
were prepared, and photographs were taken. Ampicillin was 
used as the positive control in the experiment. 


RESULTS 


cDNA cloning and characterization of amurin-9K Y 


A 279-bp cDNA sequence (GenBank accession No.: 


JX421759) encoding a 62-aa precursor was cloned from 


200 www.zoores.ac.cn 


the constructed skin cDNA library of R. kunyuensis. The 
nucleotide sequence of the cDNA and the deduced amino 
acid sequence are shown in Figure 1. To confirm the 
existence and accuracy of the cloned cDNA sequence, 
two antisense specific primers, that is, RKAMP-R1 
(5’-GCCAAGACACCCGATGTGTATTTAG-3’) and RKAMP-R2 
(5’-CCCTTTTCCACATTTTCTGGTAATT-3’), were designed 
according to the cloned cDNA sequence and coupled with the 
sense specific primer RanaAMP for PCR (synthesized skin 
cDNA library of R. kunyuensis was used as a template). The 
PCR products were sequenced and verified the identity of the 
cDNA sequence above (data not shown). 


atgttcaccttgaagaaatccctgttgcttcttttcttccttgggatcatcaacgtatct 60 
MPP Le KE SOL Be BD BO ДЖО dox eugxus 20 
ctctgtgaggaagagagaaatgccgaggaagaaagaagagatgatcccgaagaaagggcc 120 
LCEEE RN AEEE R R DDP EE RA 40 
gttgaggtggaaaaaagatttttaccattttttgcagcatgtgcaattaccagaaaatgt 180 
V EV EK R|F LP FF AACA ITREKEKC 60 
ggaaaatgatttttctaaatacacatcgggtgtcttataaaaaataaagatgaagcctac 240 
G K|* 62 
aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 219 


Figure 1 Nucleic acid sequence of cDNA encoding the 


























precursor of amurin-9K Y 
The signal peptide domain is indicated in italics, the putative mature peptide 
is boxed and the stop codon is indicated by an asterisk (*). 


The 62-aa AMP precursor deduced from the cDNA 
sequence included an N-terminal signal peptide domain of 
22 aa, an acidic amino acid residue-rich (Asp and Glu) 
Spacer peptide domain and a mature peptide domain of 16 aa 
residues. There was a conserved dibasic cleavage site Lys-Arg 
(K-R) between the spacer peptide domain and mature peptide 
domain, which was likely cleaved by trypsin-like proteases 
followed by mature peptide release. Sequence alignment of the 
mature peptide with the NCBI protein database revealed that it 
possessed highest primary structure similarity with amurin-9AM 
(GenBank accession No.: AEP84587). Hence, it was designated 
amurin-9KY according to the origin of R. kunyuensis. 

Amurin-9KY was also found to possess high sequence 
similarity to the temporin AMP family, especially the C-terminal 
Gly-Lys (G-K) di-residue, which is regarded as an amidation site 
(Figure 2). However, unlike the temporin family, amurin-9KY 
exhibited two cysteine residues forming an intra-molecular 
disulfide bridge, which appeared as a C-terminal heptapeptide 
ring (Rana box domain). The physical and chemical parameters 
of amurin-9KY were computed by ProtParam (http://web.expasy. 
org/protparam/), demonstrating a molecular weight of 1 584.97, 
theoretical pl of 8.96 and net charge of +2. 

To evaluate the influence of C-terminal amidation and the 
heptapeptide ring on the biological activities of the peptide in 
vitro, three derivatives (amurin-9KY1—3) were designed and 
chemically synthesized, with their biological functions then 
examined (Table 1). 


Antimicrobial and hemolytic assays 
As listed in Table 2, amurin-9KY was active against the tested 


gram-positive bacteria S. aureus ATCC25923, S. aureus 
090223* (IS) and N. asteroids 090312* (IS), and showed 
the strongest potency against S. aureus ATCC25923, with 
an MIC of 4.68 ug/mL. However, amurin-9KY did not show 
any antimicrobial activity against the tested gram-negative 
bacteria or fungi, consistent with the features of the temporin 
family (Conlon, 2004). The derivative amurin-9KY1 exhibited 


less potent activities against all three gram-positive bacteria. 


Interestingly, unlike amurin-9KY, amurin-9KY1 was able to kill the 
eukaryotic microbe Slime mold 090413 (IS) with an MIC value of 
75 ug/mL. However, derivatives amurin-9KY2 and amurin-9KY3 
exhibited no activities against any of the tested microbes. 

The hemolytic activities of amurin-9KY and the three 
derivatives were also examined. All four peptides showed 
slight hemolytic activity against fresh human erythrocytes at 
concentrations of 100 wg/mL, with hemolytic rates of 2%, 
15.4%, 17.9% and 20.8%, respectively. 


Table 1 Structural parameters of amurin-9KY and its derivatives 


Kunyuenin —-FLPFFAACAITRKCGK 16 К. kunyuensis 
Temporin-SHe FLPALAGIAGLLGKIFGK 18  R.saharica 
Temporin-LTd -——-FLPGLIAGIAKMLGK 15 ХА. latouchii 
Тетрогіп-1Ја —-ILPLVGNLLNDLLGK 15 КА. japonica 

Temporin-1V ——-FLPLVGKTLSGLIGK 15  R.versabilis 
Temporin-G —--FFPVIGRILNGILGK 15 R. temporaria 
Temporin-10al -—-FLPLLASLFSRLLGK 15  R.ornativentris 
Temporin-CDYb —-ILPILAPLIGGLLGK 15 R. dybowskii 











жж 
Figure 2 Amino acid sequence comparison of amurin-9KY with 


amurin-9AM from R. amurensis and temporin family AMPs 
from other ranid frogs 

All AMP sequences used were derived from the NCBI protein database. The 
two cysteine residues of amurin-9KY are indicated in gray shadow. The 
length of the AMPs is shown and species from which the AMPs are derived 


are on the right. Asterisks (*) indicate identical residues. 





Sample Amino acid sequence Molecular weight Structure characteristics 

Amurin-9KY | FLPFFAACAITRKC-NH2 1 584.96 C-terminal amidation and C-terminal heptapeptide ring C8-C14 
Amurin-9KY1 | FLPFFAACAITRKC-NH» 1 586.96 C-terminal amidation 

Amurin-9KY2 | FLPFFAACAITRKC 1 585.96 C-terminal heptapeptide ring C8-C14 

Amurin-9KY3 | FLPFFAACAITRKC 1 587.96 None 


Table 2 Antimicrobial activities of amurin-9KY and its 
derivatives 





MIC (ug/mL) 





Amurin- Amurin- Amurin- Amurin- 


Mi " 
icroorganism 9KY 9KY1 9KY2 9KY3 





Gram-positive 


Staphylococcus aureus ATCC25923 4.68 37.5 ND ND 
Staphylococcus aureus 090223* (IS) 37.5 75 ND ND 
Nocardia asteroids 090312* (IS) 37.5 75 ND ND 
Gram-negative 

Escherichia coli АТСС25922 ND ND ND ND 
Klebsiella pneumonia 1368 (IS) ND ND ND ND 
Pseudomonas aeruginosa ATCC27853 ND ND ND ND 
Fungi 

Candida albicans ATCC2002 ND ND ND ND 
Slime mold 090413 (IS) ND 75 ND ND 





MIC: Minimum inhibitory concentration. These concentrations 
represent mean values of three independent experiments performed 
in duplicate. IS: Clinically isolated strain. ND: No detectable activity 
in inhibition zone assay at a dose of 2 mg/mL. 


Anti-oxidant activity 

To date, many peptides exhibiting anti-oxidant activity have 
been identified from several species of ranid frogs, which 
constitute the excellent skin anti-radiation defense system of 
ranid frogs (He et al., 2012; Liu et al., 2010; Lu et al., 
2010; Yang et al., 2009). Most are dual-functional peptides, 
possessing both anti-oxidant and antimicrobial activities (He et 


al., 2012; Liu et al., 2010; Lu et al., 2010; Yang et al., 2009). 
In the present study, the anti-oxidant activities of amurin-9KY 
and the three derivatives were also evaluated (Figure 3). 
Amurin-9KY exhibited slight concentration-dependent DPPH 
radical scavenging activity, with an optimal S% value of 30.6% 
at a concentration of 400 ug/mL. In contrast, the derivatives 
amurin-9KY1 and amurin-9KY3 exhibited strong anti-oxidant 
activities, with S% values exceeding 60% at concentrations 
as low as 50 ug/mL. Compared to the other three peptides, 
amurin-9KY2 exhibited the lowest anti-oxidant activity, with an 
S% value of 20% at a concentration of 400 ug/mL. 


Solution structures of amurin-9KYs 


The CD spectra of amurin-9KY in water showed a negative 
band at 200 nm, indicating a random-coil conformation. In the 
membrane-mimetic solvent (50% TFE-water) the presence of 
one positive band (190 nm) and two large negative dichroic 
bands at 208 and 222 nm (-50 mdeg) in amurin-9KY was 
consistent with the o-helical conformations (Figure 4A). For 
amurin-9KY1-3, the negative dichroic bands at 208 and 222 
nm were much smaller, about —10 mdeg. The ao-helical 
structure of most active AMPs is thought to be responsible for 
the formation of pores in the membranes of target organisms, 
thus disrupting metabolic activity (Brogden, 2005). The CD 
results support the concept that amurin-9KY most likely killed 
the bacteria through membrane disruption. In addition, the 
much smaller negative 208 and 222 nm dichroic bands were 
in good agreement with the antimicrobial activity data for 
amurin-9KY 1-3. 
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Effects of amurin-9KY on microbial membrane morphology 
by SEM 

Previous studies have demonstrated that AMPs achieve 
antimicrobial activity by disrupting various key microorganism 
cell processes, with some AMPs possessing multiple 
mechanisms (Brogden, 2005; Nguyen et al., 2011). There 
are many specific targets in microbial cells for AMPs, including 
external proteins, outer surface lipids, outer membrane proteins 
(gram-negative bacteria), inner membrane, integral membrane 
proteins, nucleic acids and intracellular proteins (Liu et al., 
2017, 2018; Nguyen et al., 2011). Among them, disrupting 
the integrity of the microbial inner membrane is the most 
common mode for AMPs (Liu et al., 2017, 2018; Nguyen et al., 
2011), the disruption of which results in obvious morphological 
alteration. Here, the antimicrobial results showed that S. 
aureus ATCC25923 was most sensitive to amurin-9KY, and 
therefore it was selected to examine the induced membrane 
morphological alterations by SEM. As illustrated in Figure 5, 
the untreated S. aureus cells exhibited normal shape and 
smooth surfaces (Figure 5A). In contrast, after treatment with 
amurin-9KY for 30 min, the cellular shape and surface of S. 
aureus exhibited obvious alterations (Figure 5B). The bacterial 
cells showed rough surfaces, prevalent membrane vesicles, 
and cellular fragments, implying that amurin-9KY might act on 
the bacterial membrane and induce disruption of membrane 
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Figure 4 CD spectra of the four peptides in different solutions 
A: Amurin-9KY; B-D: Amurin-9KY1-3. 
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Amurin-9KY2 


integrity. Ampicillin-treated S. aureus exhibited no significant 
morphological alteration compared with the untreated bacteria 
(Figure 5C), except for the appearance of many granules on 
the surfaces. Ampicillin usually kills bacteria by disturbing the 
synthesis of the cell wall. 
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A: S. aureus ATCC25923 control; B: Amurin-9KY-treated S. aureus ATCC25923; C: Ampicillin-treated S. aureus ATCC25923. 


DISCUSSION 


Among the AMP families identified from ranid species, the 
C-terminal heptapeptide ring (Rana box domain, two cysteine 
residues connected by a disulfide bridge) is a common 
structural feature, which broadly exists in families such as 
brevinin-1, brevinin-2, esculentin-1, esculentin-2, ranatuerin-1, 
ranalexin, japonicin-1, nigrocin-2, odorranain-C, odorranain-D, 
odorranain-F, odorranain-G, odorranain-H, odorranain-P1 and 
hainanenin (Conlon, 2004; Duda et al, 2002; He et 
al, 2012; Li et al, 2007; Matutte et al, 2000). In 
addition, other cyclic ranid AMPs have also been identified, 
including those with a C-terminal octapeptide ring (japonicin-2), 
C-terminal hexapeptide ring (ranatuerin-2, amolopin-6), middle 
heptapeptide ring (palustrin-2) or ring formed by more than 
seven amino acids (odorranain-A, odorranain-B, odorranain-J, 
odorranain-P2, odorranain-T, odorranain-U and ranacyclins) 
(Conlon, 2004; Duda et al., 2002; He et al., 2012; Li et al., 2007; 
Matutte et al., 2000). To date, three AMP families characterized 
with C-terminal amidation have been identified from ranid frogs, 
including temporins, ranacyclins and tiannanenin (Conlon, 
2004; He et al., 2012; Li et al., 2007). Among the AMP families 
identified from ranid frogs so far, ranacyclin is a unique family 
possessing both a disulfide bridge-formed ring and amidated 
C-terminus (Mangoni et al., 2003). The ring structure of 
ranacyclins is composed of 11 residues and is located in the 
middle of the sequences. Comparatively, the novel AMP in the 
present study (amurin-9KY) is the first reported to have both a 
C-terminal heptapeptide ring and amidated C-terminus. 


The three derivatives (amurin-9KY 1-3) were designed here 
to evaluate the influence of C-terminal amidation and the 
C-terminal heptapeptide ring on the biological activities of 
amurin-9KY. The four peptides were chemically synthesized 
and their in vitro functions were subsequently examined. 
The antimicrobial assay results indicated that C-terminal 
amidation played an important role in the antimicrobial activity 
of amurin-9KY, whereas the heptapeptide ring contributed 
no obvious influence. SEM demonstrated that amurin-9KY 
induced obvious bacterial membrane morphological alteration, 
indicating that it might act through the disruption of bacterial 
membrane integrity. 


Amurin-9KY possessed strong antimicrobial activity against 
gram-positive bacteria and low hemolytic activity, consistent 


with the features of temporin family AMPs. Previous study 
of temporin-1Od from the Japanese mountain brown frog R. 
ornativentris demonstrated that it possessed high antimicrobial 
potency towards S. aureus due to a positive charge associated 
with the free N-terminal amino group (Kim et al., 2001). The 
current study further confirmed that the free amino group is 
essential for the antimicrobial activity of AMPs, no matter at 
which terminus (N or C) it is located. Previous structure-activity 
analysis of brevinin 1E, a brevinin-1 family AMP identified 
from R. esculenta, demonstrated that the elimination of the 
intra-disulfide bridge did not greatly affect the antimicrobial 
activity (Kwon et al., 1998), which was further verified in the 
present result. 


Regarding hemolytic activity, both C-terminal amidation 
and the heptapeptide ring significantly reduced the hemolytic 
activity of amurin-9KY compared with the derivatives without 
these two structural features. Elimination of the C-terminal 
heptapeptide ring greatly increased the anti-oxidant activity 
of amurin-9KY, whereas elimination of C-terminal amidation 
did not affect it at all. These results are likely because the 
unpaired thiol group of cysteine generated after elimination 
of the heptapeptide ring acted as an electron donator to the 
radical, which is crucial for the anti-oxidant activity of peptides 
(Akerstróm et al., 2007; Liu et al., 2010). 
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ABSTRACT 


Protease inhibitors have been reported rarely from 
the leech Hirudinaria manillensis. In this study, we 
purified a novel protease inhibitor (bdellin-HM-2) with 
anticoagulant properties from H. manillensis. With a 
molecular weight of 1.4x10^, bdellin-HM-2 was also 
characterized with three intra-molecular disulfide 
bridges at the N-terminus and multiple HHXDD and 
HXDD motifs at the C-terminus. cDNA cloning revealed 
that the putative nucleotide-encoding protein of 
bdellin-HM-2 contained 132 amino acids and was 
encoded by a 399 bp open reading frame (ORF). 
Sequence alignment showed that bdellin-HM-2 shared 
similarity with the "non-classical" Kazal-type serine 
protease inhibitors, but had no inhibitory effect on 
trypsin, elastase, chymotrypsin, kallikrein, factor Xlla 
(FXlla), factor Xla (FXla), factor Xa (FXa), thrombin, 
or plasmin. Bdellin-HM-2 showed anticoagulant effects 
by prolonging the activated partial thromboplastin time 
(aPTT), indicating a role in enabling H. manillensis to 
obtain a blood meal from its host. Our results 
suggest that bdellin-HM-2 may play a crucial role in 
blood-sucking in this leech species and may be a 
potential candidate for the development of clinical 
anti-thrombotic drugs. 


Keywords: Hirudinaria manillensis; Bdellin-HM-2; 
"Non-classical" Kazal inhibitors; Blood sucking; 
Anticoagulant; Anti-thrombotic drugs 


INTRODUCTION 
Protease inhibitors occur naturally in living organisms, including 


animals (Shadrin et al., 2015; Vicuna et al., 2015; Wang et al., 
2005; Zhang, 2006), plants (Kim et al., 2009; Ryan, 1990), 


Science Press 


fungi (Sabotic & Kos, 2012), and bacteria (Supuran et al., 
2002). They have multifunctional roles in many physiological 
processes and play an important role in biological functions of 
venomous animals, such as in predation (Birrell et al., 2007) 
and defense (Ali et al., 2002). To prevent clotting during blood 
feeding from a host, hematophagous animals have developed 
various mechanisms to interfere with blood coagulation 
(Markwardt, 1996). Among the inhibitors involved in coagulation, 
protease inhibitors are the most prominent anticoagulants 
currently described and characterized from leeches (Dodt, 
1995) and insects (De Marco et al., 2010; Mende et al., 1999). 

There are at least four types of protease inhibitors, including 
serine, cysteine, aspartic, and metalloprotease inhibitors (Leung 
et al., 2000). The Kazal family is one of the best-known 
groups of serine protease inhibitors (Rimphanitchayakit & 
Tassanakajon, 2010). Kazal-type inhibitory can be sorted into 
classical and non-classical Kazal domains. The classical 
Kazal domain has two residues between cys4 and cys5, 
whereas the non-classical Kazal inhibitor has a spacer region 
between cys4 and cys5, ranging from three to seven residues 
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(Moser et al., 1998). There are highly homologous three- 
dimensional structures in the Kazal-type serine proteinase 
inhibitors regardless of length of amino acid sequences 
between the cysteines and amino acid sequence variation 
(Eigenbrot et al., 2012). The P1 residue, located in the second 
amino acid downstream of the second conserved cysteine 
residue, is inserted into the S1 specificity pocket of the 
protease in a substrate-like way (Bode & Huber, 1992; 
Laskowski & Kato, 1980). 

Several Kazal-type serine protease inhibitors have been 
characterized from leeches. A few “non-classical” Kazal 
inhibitors have been reported from different leeches, including 
bdellin-B-3 (Fink et al., 1986), bdellin-KL (Kim et al., 2001), 
and bdellin-HM (Lai et al., 2016). In this study, bdellin-HM-2 
was purified and characterized from the leech H. manillensis. 
To the best of our knowledge, bdellin-HM-2 is the first Kazal- 
type serine protease inhibitor displaying anticoagulant 
properties identified from H. manillensis. 


MATERIALS AND METHODS 


Collection of crude extracts 

The H. manillensis leeches were purchased from Jinbian 
aquafarm, Qinzhou City, Guangxi Province in China. The 
leeches were still alive when transported to the laboratory. We 
prepared the crude extracts from the leech heads as 
described previously (Lai et al., 2016). In short, leech heads 
were separated from the bodies, washed in 0.9% saline, 
quickly frozen, and then ground in liquid nitrogen. 


Purification of bdellin-HM-2 

Purification of bdellin-HM2 followed similar methods described 
in our previous published article (Lai et al., 2016). Briefly, 
crude extracts were dissolved in 50 mmol/L Tris-HCI buffer 
(pH 8.9) and subsequently separated by a DEAE Sephadex A- 
50 column (5 cm diameter, 60 cm length, GE, USA). Elution 
was performed at a flow rate of 15 mL/h at 4 °C and 3.0 mL 
fractions were collected in separate tubes. The absorbance of 
the fractions was monitored at both 215 and 280 nm. Fractions 
that could prolong the activated partial thromboplastin time 
(aPTT) were pooled and lyophilized prior to further purification. 
The powder from the previous step was dissolved and loaded 
for reverse-phase high-performance liquid chromatography 
(RP-HPLC) on a C,, column (Waters, Milford, MA, USA, 5 um 
particle size, 250 mmx4.6 mm). Elution was carried out with a 
linear gradient of 1096 — 60% solution B (99.9% acetonitrile, 
0.196 TFA) for 60 min at a flow rate of 1 mL/min. The eluted 
fraction that prolonged aPTT was collected. 


Mass spectrometric analysis and peptide sequencing 

The molecular weight of the collected fraction was analyzed 
by matrix-assisted laser desorption ionization time-of-flight 
mass spectrometry (MALDI-TOF-MS, Autoflex speed TOF/ 
TOF, Bruker Daltonik GmbH, Bruker Corporation, Germany) 
using positive ion and linear mode, with specific operating 
parameters including a 20 kV ion acceleration voltage, 50- 
time accumulation for single scanning, and 0.196 accuracy of 
mass determinations. The partial peptide sequence of the N- 
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terminal was determined by automatic Edman degradation on 
a pulsed liquid-phase sequencer (PPSQ-31A, Shimadzu 
Corporation, Japan). 


RNA extraction and cDNA library construction 

Total RNA from the head of H. manillensis was extracted using 
Trizol reagent (Life Technologies, Carlsbad, CA, USA) according 
to the manufacturer' s instructions and was dissolved in 
RNase-Free water. A SMART ™ PCR cDNA construction kit 
(Clontech, Palo Alto, CA, USA) was used for synthesizing 
CDNA, as described previously (Lai et al., 2016). 


Screening of cDNA encoding bdellin-HM-2 

To screen the cDNA encoding the precursor of bdellin-HM-2, 
synthesized cDNA was used as the template for PCR, 
following previously described methods (Lai et al., 2016). 
Briefly, two pairs of oligonucleotide primers (Table 1) were 
used in PCR reactions, where primers 1 and 3 were designed 
according to the partial N-terminal sequence of bdellin-HM-2 
and primers 2 and 4 were from the SMART ™ PCR cDNA 
construction kit. The PCR conditions were as described 
previously (Lai et al., 2016). 


Table1 Primers used for cDNA cloning of bdellin-HM-2 





Primer Sequence (5'-3') 

1 AACAGGTTTGCGGAAGT 

2 AAGCAGTGGTATCAACGCAGAGT 
3 AATTCCAGGGTACAGACG 

4 ATTCTAGAGGCCGAGGCGGCCGA 


Primers 1 and 2 for signal peptide cloning; primers 3 and 4 for mature 
peptide cloning. 


Effects of bdellin-HM-2 on blood coagulation 

For aPTT assay, the aPTT reagent (50 uL, F008-1, Nanjing 
Jiancheng Bioengineering Institute, China) was incubated with 
50 uL of plasma and different concentrations of bdellin-HM-2 
(0.7 and 1.4 pmol/L). After 3-min incubation, CaCl, (50 uL, 
25 mmol/L) preheated at 37 °C for 5 min was added, and the 
clotting curve was monitored at 650 nm using an enzyme- 
labeled instrument (Epoch BioTek, USA) for 2 min. To test the 
prothrombin time (PT), plasma (50 uL) was incubated with 
different concentrations of bdellin-HM-2 (0.7 and 1.4 pmol/L) 
for 3 min at 37 °C, after which the PT reagent (100 uL, F007, 
Nanjing Jiancheng Bioengineering Institute, China) preheated 
at 37 °C for 15 min was added and the clotting curve was 
monitored at 650 nm using the enzyme-labeled instrument for 
30 s. 


Effects of bdellin-HM-2 on proteases 

Effects of bdellin-HM-2 on proteases, including trypsin, elastase, 
chymotrypsin, kallikrein, factor XIla (FXlla), factor Xla (FXla), 
factor Xa (FXa), thrombin, and plasmin were tested using the 
corresponding chromogenic substrates. The testing enzyme 
was incubated with different concentrations (0, 0.7, 1.4, 2.8, 
and 5.6 umol/L) of bdellin-HM-2 in 60 uL of 50 mmol/L Tris 
buffer (pH 7.4) for 5 min, with a certain concentration of 


chromogenic substrate then added. Absorbance at 405 nm 
was monitored immediately and the kinetic curve was recorded 
using the enzyme-labeled instrument for 30 min. Bovine 
pancreas trypsin, elastase, chymotrypsin, and plasmin were 
all obtained from Sigma (USA) and the enzyme concentrations 
used were 800, 400, 400, and 20 nmol/L, respectively. The 
corresponding chromogenic substrates (Sigma, USA) were 
Gly-Arg-p-nitroanilide dihydrochloride for trypsin, N- 
methoxysuccinyl-Ala-Ala-Pro-Val-p-nitroanilide for elastase, N- 
succinyl-Gly-Gly-Phe-p-nitroanilide for chymotrypsin, and Gly- 
Arg-p-nitroanilide dihydrochloride for plasmin. The 
concentration of all substrates in the reactions was 0.2 mmol/L. 
The concentrations used for kallikrein, FXla, and FXa (Enzyme 
Research Laboratory, USA) were 400, 400, and 20 nmol/L, 
respectively, and the corresponding chromogenic substrates 
were H-D-Pro-Phe-Arg-pNA : 2HCI (Hyphen Biomed, France), 
H-D-Pro-Phe-Arg-pNA : 2HCI (Hyphen Biomed, France), and 
CH,0CO-D-CHA-Gly-Arg-pNA-AcOH (Sigma, USA), 
respectively. The concentration of all three substrates in the 
reaction was 0.2 mmol/L. Human thrombin (Sigma, USA, 
10 nmol/L) and FXIla (Enzyme Research Laboratories, USA, 
10 nmol/L) were reacted with 0.2 mmol/L chromogenic 
substrate of H-D-Phe-Pip-Arg-pNa · 2HCI (Hyphen Biomed, 
France) and H-D-Pro-Phe-Arg-pNA : 2HCI (Hyphen Biomed, 
France), respectively. 


RESULTS 


Purification of bdellin-HM-2 

The crude extracts from H. manillensis were resolved into 
several fractions by DEAE Sephadex A-50 column. The 
fraction that prolonged the aPTT was indicated by a bar 
(Figure 1A). We then obtained the purified peptide exerting an 
aPTT inhibitory effect, named bdellin-HM-2 (indicated by an 
arrow in Figure 1B), using a С, RP-HPLC column. MALDI- 
TOF-MS showed that bdellin-HM-2 had a molecular weight 
(MW) of 14141.5 (Figure 1C). 


Primary structure of bdellin-HM-2 

Based on automatic Edman degradation, the partial N- 
terminal sequence of bdellin-HM-2 was determined to be 
ETECVCTLELKQVCGS. According to the N-terminal sequence, 
degenerate primers were designed (Table 1) to clone the 
cDNA encoding the precursor of bdellin-HM-2 from the cDNA 
library. A 399 bp cDNA encoding the precursor of bdellin-HM-2 
was obtained. The cDNA had an open reading frame (ORF) of 
396 nucleotides coding a pro-protein of 132 amino acids, 
including a signal peptide of 18 residues (indicated by box) 
and mature bdellin-HM-2 of 114 residues (Figure 2A). The 
theoretical MW of mature bdellin-HM-2 was 13144.78, which 
was not consistent with the observed mass by mass 
spectrometry analysis (Figure 1C). This inconsistency may be 
due to post-translational modification of the protein. Sequence 
alignment showed similarity to bdellin-KL (Kim et al., 2001), 
bdellin-B-3 (Fink et al., 1986), and bdellin-HM (Lai et al., 
2016), which are ‘non-classical’ Kazal serine protease inhibitors 
(Figure 2B). Multiple sequence alignment showed that the six 
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Figure 1 Purification of bdellin-HM-2 from H. manillensis 

A: Crude extracts were fractionated using DEAE Sephadex A-50 ion 
exchange. Fraction exerting inhibitory activity on aPTT is indicated by 
a bar (—). B: Fraction exerting inhibitory activity on aPTT was further 
purified by C,, reverse-phase high-performance liquid chromatography 
(RP-HPLC). Protein peak exerting inhibitory activity on the aPTT is 
marked by an arrow. C: Matrix-assisted laser desorption ionization 
time-of-flight (MALDI-TOF) analysis of purified native bdellin-HM-2. 


cysteine residues and threonine-tyrosine residues were highly 
conserved among different species (Figure 2C). There were 
multiple HHXDD and HXDD motifs at the C-terminus of bdellin- 
HM-2. 


Anticoagulant activity of bdellin-HM-2 

Underthe assay conditions, bdellin-HM-2 exerted anticoagulatory 
activity by inhibiting aPTT (Figure 3A), whereas no inhibitory 
activity was observed on PT (Figure 3B). Compared with the 
control with an aPTT of ~60 s, the aPTT was prolonged to ~ 
100 s after 0.7 and 1.4 pmol/L bdellin-HM-2 treatment, 
suggesting that bdellin-HM-2 acts on the intrinsic pathway. 
Bdellin-HM-2 had no effect on trypsin, elastase, chymotrypsin, 
kallikrein, FXlla, FXla, FXa, thrombin, or plasmin (Figure ЗС). 
All enzyme activity test results were plotted (Figure 4). 
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atgaagtttttgatcgctcttctctttctctgctctctggtcagcatcgaggctgagaca 

MK FLITALLFLCSLVSTEAET 20 
gaatgcgtctgtaccctggaattgaaacaggtttgcggaagtgatggtgtcacatatgac 
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Figure 2 cDNA sequence encoding bdellin-HM-2 precursor and 
sequence alignment with other protease inhibitors 

A: Nucleotide sequence encoding bdellin-HM-2 precursor and deduced 
amino acid sequence. Signal peptide is boxed. Bar (-) indicates stop 
codon. B: Sequence comparison of bdellin-HM-2 with bdellin-HM, 
bdellin-KL, and bdellin-B-3. Identical amino acid residues are indicated 
by an asterisk (*). C: Multiple sequence alignment of Kazal domain 
from H. manillensis (bdellin-HM-2), Aedes aegypti (AaKPI ABF18209), 
Anemonia sulcata (AsEI 1Y1B), and H. medicinalis (СОТ! P80424). 
Conserved threonine-tyrosine residues between cysteine 3 and 4 are 
indicated. Conserved cysteine motifs are also indicated. 


DISCUSSION 


Several protease inhibitors exerting anticoagulant effects have 
been found from leeches (Hong & Kang, 1999; Markwardt, 
2002; Salzet et al., 2000; Salzet, 2001; Strube et al., 1993; 
Tuszynski et al., 1987). In this report, a novel protease 
inhibitor (bdellin-HM-2) with anticoagulant effects was purified 
and further characterized from H. manillensis for the first time. 
The cDNA encoding bdellin-HM-2 precursor was cloned from 
the cDNA library, and the mature bdellin-HM-2 consisted of 
114 amino acid residues. MALDI-TOF-MS showed that the 
MW of bdellin-HM-2 was 14141.5, compared to the theoretical 
molecular weight of 13144.78, a difference of 996.72, which is 
not consistent with the theoretical value. Research shows 
glycosylation influences the function of protein, governs 
physiology, and contributes to disease (Ohtsubo & Marth, 
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Figure 3 Effects of bdellin-HM-2 on aPTT 
Bdellin-HM-2 action on aPTT (A) and PT (B). C: Bdellin-HM-2 effects 
on trypsin, elastase, chymotrypsin, kallikrein, FXlla, FXla, FXa, 
thrombin, and plasmin were analyzed by hydrolysis of chromogenic 
substrates. Data are means+SD of six independent experiments. 
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2006). We speculated that bdellin-HM-2 was O-glycosylated 
at Thr-20, Thr-25, Ser-34, Thr-38, and Thr-46. (Gupta & Brunak, 
2002), although further research on these O-glycosylation 
sites should be performed in the future. 

Kazal-type inhibitors with one or more Kazal domains are 
characterized by multiple HHXDD and HXDD motifs in their 
amino acid sequences (Laskowski & Kato, 1980) and by their 
highly homologous three-dimensional structures (Van et al., 
1995). Each Kazal domain usually contains six conserved 
cysteine residues forming three intra-molecular disulfide 
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Figure 4 Bdellin-HM-2 had no effect on proteases 


Bdellin-HM-2 effects on trypsin (A), elastase (B), chymotrypsin (C), kallikrein (D), FXlla (E), FXla (F), FXa (G), thrombin (H), and plasmin (1). Data 
represent at least six independent experiments and are presented as means+SD. 


bridges (Laskowski & Kato, 1980; Magert et al., 1999). P1 
residue, which contributes to the inhibitory specificity, is 
located at the second position after the second cysteine 
residue of the Kazal domain (Bode & Huber, 1992). Although 
bdellin-HM-2 showed high similarity to bdellin-HM and bdellin- 
KL by sequence analysis and belongs to the family of non- 
classical Kazal domains, enzyme activity tests showed that 
bdellin-HM-2 had no inhibitory effects on trypsin, elastase, 
chymotrypsin, kallikrein, FXlla, FXla, FXa, thrombin, or plasmin 
under the assay conditions. Sequence alignment showed that 
the P1 residue of bdellin-HM-2 was different from bdellin-HM, 
bdellin-KL, and bdellin-B-3. The difference in P1 residue was 
considered the cause of the enzyme activity test results. 
Bdellin-HM-2 prolonged the aPTT, implying that bdellin-HM- 
2 functioned to help H. manillensis obtain a blood meal by 
inhibiting blood coagulation. Results showed that the activity 
was dose-independent. Further work to identify the target of 
bdellin-HM-2 in blood is necessary. Blood-sucking animals 
obtain a blood meal by overcoming the host’ s natural blood 
coagulation (De Marco et al., 2010; Dodt, 1995; Markwardt, 
1996; Mende et al., 1999). The anticoagulant peptide obtained 
from H. manillensis not only facilitates our understanding of 
the mechanism of blood feeding for H. manillensis, but also 
provides a new candidate for the development of clinical 


anticoagulant drugs. 

In conclusion, bdellin-HM-2 identified from H. manillensis 
prolonged the aPTT but exhibited no influence on PT and no 
inhibitory activity on trypsin, elastase, chymotrypsin, kallikrein, 
FXlla, FXla, FXa, thrombin, or plasmin under the assay 
conditions. Further research on O-glycosylation sites will be 
performed in the future. Bdellin-HM-2 is the first identified 
Kazal-type serine protease inhibitor from H. manillensis that 
shows a potent anticoagulant effect. 
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ABSTRACT 


Accurate information on eye position in the orbit 
is available from visual feedback, efference copy 
of the oculomotor commands and proprioceptive 
signals from the extraocular muscles (EOM). Whereas 
visual feedback and oculomotor commands have 
been extensively studied, central processing of EOM 
proprioceptive signals remains to be elucidated. A 
challenge to the field is to develop an approach to 
induce passive eye movements without physically 
contacting the eyes. A novel method was developed 
to generate passive eye movements in rats. A small 
rare-earth magnet disk (0.7 mm diameter, 0.5 mm 
thickness) was attached to the surface of a rat's 
eyeball. A metal rod (5 mm diameter) wrapped 
with an electromagnetic (EM) coil was placed near 
the magnet (8-15 mm). By passing currents to the 
EM coil, electromagnetic force (EMF) was generated 
and acted upon the magnet and induced passive 
eye movements. The EMF induced well-defined 
passive eye movements, whose directions were 
dependent on current polarity and amplitudes and 
peak velocities were dependent on current intensity 
and duration. Peak velocities of the EMF-induced 
eye movements were linearly related to amplitudes, 
exhibiting main sequence relationships similar to that 
of saccades in awake rats and eye movements induced 
by electrical microstimulation of the abducens nucleus 
in anesthetized rats. Histological examination showed 
that repetitive EMF stimulations did not appear to result 
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in damages in the EOM fibers. These results validated 
the EMF approach as a novel tool to investigate 
EOM proprioceptive signals and their roles in visual 
localization and gaze control. 


Keywords: Eye movement; Proprioception; Extraocular 
muscles; Stretch reflex 


INTRODUCTION 


An important issue in modeling visual localization and gaze 
control is to determine whether the central nervous system 
(CNS) obtains eye position information in the orbit via efference 
copy (or corollary discharge) of ocular motor commands or 
proprioception signals from the extraocular muscles (EOM) 
or the orbital connective tissues. Whereas Helmholtz (1925) 
emphasized the roles of efference copy, Sherrington (1918) 
emphasized the roles of EOM proprioception. The two camps 
have been debating this issue over the past century. On one 
hand, earlier studies showed that humans are aware of the 
passive displacements of the eyes in darkness (Skavenski, 
1972) and altering this proprioceptive information consistently 
leads to errors in visual localization (Allin et al., 1996; Balslev 
& Miall, 2008; Bridgeman & Stark, 1991; Campos, 1986; 
Campos et al., 1989; Gauthier et al., 1990; Lennerstrand 
et al., 1997). Recent finding of EOM proprioceptive signals 
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in the somatosensory cortex (Wang et al., 2007) and its 
roles in computing spatial maps (Xu et al., 2012) provided 
compelling evidence that the CNS utilizes the proprioceptive 
eye position signals. On the other hand, bilateral section 
of the monkey trigeminal nerves did not reduce accuracy of 
open-loop pointing (Lewis et al., 1994) and memory-guided 
saccades (Guthrie et al., 1983). Based on the assumption that 
the trigeminal nerves carry eye proprioceptive information to 
the brain (Porter et al., 1983), these results suggested that eye 
proprioception signal is not necessary for visual localization 
and gaze control because efference copy of the oculomotor 
motor commands is sufficient to implement these tasks. 

The lack of understanding of the EOM proprioceptive signals 
is partially attributable to the fact that location of the peripheral 
sensors and their connections to the brain are still a matter 
of debate (Billig et al., 1997; Lienbacher et al., 2011a, 201 1b; 
Zimmermann et al., 2011, 2013). However, it is also attributable 
to limitations of the three types of methods that have been 
used to induce passive eye movements. One method is to 
eliminate the EOM proprioceptive signals by use of surgical 
(Guthrie et al., 1983; Lewis et al., 1998) or pharmacological 
approaches (O’Keefe & Berkley, 1991). The second method 
is to mechanically displace an eye with a suction lens (Keller 
& Robinson, 1971) or by having the subject manually press 
on one of their eyes. The third method is to perturb 
the proprioceptive representation by vibrating the extraocular 
muscle tendon or electrically stimulating extraocular muscles 
(Han & Lennerstrand, 1999). These approaches suffer from 
the fact they may disrupt ongoing behaviors and create other 
sensory or motor signals. In the present study, we described 
a novel approach, which employs electromagnetic force (EMF) 
to induce passive eye movements in rats. The goal of the 
experiments was to demonstrate that the EMF approach is an 
effective tool for producing eye movements without inducing motor 
or non-proprioceptive somatosensory signals. Consequently, it 
can be used to determine the roles of eye position proprioceptive 
signals in visual localization and gaze control. 


MATERIALS AND METHODS 


Animals 

A total of nine pigmented, female Long-Evans rats (Harlan 
Labs, Indianapolis, IN, USA) weighing 250-300 g were 
used in this study. Six of them were for EMF-induced eye 
movement experiments and three of them were for histological 
experiments. All procedures were carried out in accordance with 
NIH guidelines and approved by the Institutional Animal Care and 
Use Committee at the University of Mississippi Medical Center. 


Surgical procedures 

All surgical procedures were performed aseptically, as 
described before (Zhu et al., 2011). Briefly, a rat was implanted 
with a small head holder on the skull and allowed 7 days for 
recovery before eye movement tests. On the day of testing, 


the rat was first anesthetized with 5% isoflurane-O» (2 L/min). 


Once the rat was sedated, the isoflurane level was maintained 
at 2.5% throughout the rest of the procedure. A drop of 
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ophthalmic solution (Proparacaine Hydrochloride Ophthalmic 
Solution, USP 0.5%, Valeant Pharmaceuticals North America 
LLC, NJ, USA) was applied on the top of corneal surface for 
local anesthesia. Next, a small magnetic disk ~0.7 mm in 
diameter and 0.5 mm in thickness was attached by 3 mol/L 
Vetbond (n-butyl-cyanoacrylate) tissue adhesive (St Paul, MN, 
USA) to the surface of the eyeball nasal to the pupil (Figure 1A). 
The line connecting the magnet poles was tangential to eye 
surface. A small soft iron core rod (5 mm diameter) was 
placed within 8—15 mm of the magnet and aligned with the line 
connecting the magnet poles. The end of the rod away from 
the magnet was inside of an electromagnetic coil (Figures 1, 2), 
which generated magnetic field to magnetize the rod when 
current passed through the coil. The magnet field amplitude 
and onset were monitored by recording the EM coil current 
signal. The larger the activation current intensity (0.1—0.6 A), 
the stronger the electromagnetic field it generates. Based 
on the specifications of the EM coil (Figure 2), we used the 
following equation to estimate the magnetic flux (B) generated 
by activation current. 


В = n x 1/(Leore + Lgap/Ho) (1) 


Where B is magnetic flux in Teslas, n is number of turns, | is 
activation current in Amperes, Lcore is the length of the core 
in meters, сар is the distance from the core, u is permeability 
of soft iron and Up is permeability of air. At 10 mm from the 
core, an activation current of 0.6 A generated a magnetic flux 
of about 0.2 Teslas. 


Eye movement recording 

Horizontal and vertical eye position signals were recorded 
using a video-based ISCAN ETS-200 eye tracking system 
(ISCAN, Burlington, MA, USA). An infrared camera equipped 
with a zoom lens (Computar Optics Group, Japan) was 
attached to the platform of the rotator/sled and was focused 
on the left eye of the rat, which was secured to platform via 
the head holder. The rat's eye was illuminated by a standard 
ISCAN multiple infrared LED illuminator attached to the camera 
to produce a reference of corneal reflection (CR) for measuring 
eye position in the orbit. The eye tracker tracks the pupil center 
and the CR at a speed of 240 frames per second with a spatial 
resolution of 0.1 degrees. The differences of the two signals 
provided real-time signals related to eye position. Calibration 
was achieved by rotating the camera from left 10 degrees to 
right 10 degrees around the vertical axis of the recorded eye 
(de Jeu & de Zeeuw, 2012). 


Electrical stimulation of abducens nucleus 

Rats were anesthetized by isoflurane with their heads fixed on 
the stereotaxic frame through the use of surgically implanted 
head holder. The right occipital bone was opened and the 
cerebellum exposed. A microelectrode pulled from a thin glass 
pipette (OD: 1.2 mm; ID: 0.9 mm; Sutter Instruments, Novato, 
CA, USA) was filled with 3 mol/L sodium chloride (10—20 MQ) 
and was advanced to the abducens nucleus using coordinates 
in the atlas of Paxinos & Watson (2009). The abducens nucleus 
was identified as the site where abducting eye movements 


were evoked by low intensity electrical stimulation (10 mA, 
200 ms duration, 100 Hz, pulse duration 0.2 ms) delivered 
through the recording electrode. Stimulation-evoked eye 
movement was recorded by the eye tracker. By varying current 
intensity, electrical stimulation evoked eye movements with 
varying amplitudes were produced. After completion of the 
experiments, the recording sites were marked by ejecting fast 
green with a negative 7 mA current for 10 min, and their 
locations were subsequently verified histologically. 
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Figure 1 Schematic illustration of the experimental setup 

After anesthesia, a small magnet was attached to the surface of the rat's left 
eyeball nasal to the pupil (A). The line connecting the poles of the magnet 
was in the direction tangent to the eye surface. A small metal rod (B) was 
placed within 8—15 mm of the magnet and aligned with the line connecting 
its poles. The end of the metal rod away from the magnet was surrounded 
by an electromagnetic coil (C), which generated magnetic field to magnetize 
the metal rod when current passed through the coil. 
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Figure 2 Schematic illustration of the electromagnetic device 
It has a 59 mm soft iron core rod wrapped with ~2 600 turns of enamel 
coated magnet wire (0.25 mm diameter). The rod is attached to an aluminum 
heatsink and an aluminum support rod. 


Aluminum support rod 


Data acquisition and analysis 

Signals related to horizontal and vertical eye position and 
detecting coil voltage were sampled at 1 kHz at 16 bits 
resolution by a CED Power 1401 system (Cambridge 
Electronics Devices, Cambridge, UK). Eye movement 
responses were analyzed using Spike2 (Cambridge Electronics 
Devices, Cambridge, UK). Raw eye position data were filtered 
and differentiated with a band-pass of DC to 50 Hz to obtain 
eye velocity data. Trials in which there was no pupil tracking 
within 50 ms of the onset of the EMF onset were rejected 
(less than 5%). Trials in the data stream were aligned with the 
EMF onset and averaged (~100 trials per condition) to obtain 
low-noise estimates of eye position and velocity as a function 
of time. The EMF induced a displacement of eye position, 
whose direction was dependent on the EMF polarity (attractive 
or repulsive). After the EMF was turned off, the eye returned to 


its initial position with an exponential time course, which was 
fitted by a sum of two exponentials (Anderson et al., 2009). 


A(t) = A e Т!) + Aəe (/T2) (2) 


Examination of extraocular muscle morphology 

In order to determine whether repetitive EMF-induced 
movements of the eye damaged the orbital contents, the 
eyes of three rats were examined histologically. Following 
completion of >3 000 EMF-induced eye movements toward the 
ear (i.e., stretching the medial rectus), rats were sacrificed 
via intracardiac perfusion with 4% paraformaldehyde in pH7.2, 
0.1 mol/L phosphate-buffered saline (PBS). The orbit and 
its contents (bone, globe, extraocular muscles and glands) 
were harvested. The globe and extraocular muscles were 
placed in 30% sucrose in pH 7.2, 0.1 mol/L phosphate buffer 
(PB) solution. They were then cryoprotected by moving 
them through increasingly concentrated solutions of Cryomatrix 
embedding medium (Shandon Cryomatrix, Thermo Fisher 
Scientific Inc., DE, USA). The globe was then oriented in a 
disposable mold containing Cryomatrix and frozen. Once the 
specimen was cooled to —19 °C, individual sections containing 
the globe and extraocular muscles were cut serially at 35 um 
on a cryostat (Thermo Shandon Cryotome Е cryostat, Thermo 
Fisher Scientific Inc., DE, USA). Sections were secured to cold, 
gelatin coated glass slides. They were then stained using 
Masson Trichrome stain, cleared and cover slipped. Digital 
photographs were taken using NIS-Elements AR and a Nikon 
Eclipse E600 light microscope equipped with a 1.5 megapixel 
Nikon DS-Ri1 high resolution camera. When necessary, 
images were adjusted for brightness, contrast and color using 
Adobe Photoshop CS5 to replicate the image as it appeared 
when visualized under the microscope. 


RESULTS 


Passive eye movement induced by EMF 

The EMF induced well-defined eye movements in both 
horizontal and vertical directions. In the example shown in 
Figure 3, negative current pulses (20 ms, —0.55 A) generated 
attractive EMF to the magnet at a distance of 10 mm, which 
rotated the left eye toward the left ear with a latency of ~5 
ms, resulting in a passive eye movement of 11.4 deg in the 
leftward horizontal direction and 2.2 deg in the upward direction 
(Figure 3A). About 11 ms after the EMF was turned off, the 
elasticity of the eye plant drove the eye toward its initial position 
with exponential a time course. Plots of this movement could 
be fitted with two time constants (8.5 ms and 557 ms, R?-0.99, 
Figure 3A, blue dotted line). Positive current pulses (20 ms, 
0.55 A) generated repulsive EMF to the magnet, which rotated 
the left eye away from the left ear with a latency of ~5 ms, 
resulting in a passive eye movement of 7.4 deg in the rightward 
horizontal direction and 1.0 deg in the downward direction 
(Figure 3B). About 9 ms after the EMF was turned off, the 
elasticity of the eye plant drove the eye back to its initial position 
with an exponential time course, which could be fitted with two 
time constants (6.7 ms and 228.1 ms, R?-0.99, Figure 3B, blue 
dotted line). 
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Figure 3 Averaged eye position and velocity responses to EMF 
(Rat1) 

A: Attractive EMF induced passive eye movements of the left eye toward the 
left ear. B: Repulsive EMF induced passive eye movements of the left eye 
away from the left ear. The vertical dashed lines indicate EMF onset (current 
pulse: 0.55 A; 20 ms). Horizontal dashed lines indicate initial eye position 
and baseline eye velocity (zero). Upward traces indicate eye movements 
away from the left ear. Downward traces indicate eye movements toward the 
left ear. 


Time from EMF onset (ms) 


To assess the effectiveness of using EMF to induce passive 
eye movements, we examined effects of EM coil current 
intensity (0.1 A—0.6 A) on amplitudes of EMF-induced eye 
movement. Figure 4A shows representative eye responses 
(Rat2) at two EM coil current intensities (0.3 A and 0.4 A) and 
two polarities (repulsion and retraction) at a distance of 10 mm. 
EM coil currents with higher intensities induced eye movements 
with larger amplitudes. Figure 3B summarizes eye movement 
responses from four rats with attraction EMF (black symbols 
and lines) and three rats with repulsion EMF (grey symbols and 
lines). Eye movement amplitudes increased linearly with EM 
coil intensity. For the four cases with attraction EMF (Figure 4B, 
black lines), the slopes of eye movement-current intensity are 
23.6+0.9 deg/A, 24.3+3.6 deg/A, 10.2+1.1 deg/A and 3.9+1.4 
deg/A, respectively. For the three cases with repulsion EMF 
(Figure 4B, grey lines), the slopes of eye movement-current 
intensity are 12.2+0.9 deg/A, 9.4+0.8 deg/A and 3.2+0.2 deg/A, 
respectively. 

We also examined effects of EM coil current duration 
(5-30 ms) on amplitudes of the EMF-induced eye movement. 
Figure 5A shows representative eye responses at three EM coil 
current durations (15 ms, 20 ms and 25 ms, 0.55 A intensity) 
and two polarities (repulsion and retraction) at a distance of 
10 mm. EM coil currents with longer durations induced eye 
movements with larger amplitudes. Figure 5B summarizes 
eye movement responses from 2 rats with attraction EMF 
(black symbols and lines) and repulsion EMF (grey symbols 
and lines). Eye movement amplitude increased linearly with 
EM coil current duration. For attraction EMF, the slopes of 
eye movement-current duration are 0.59+0.03 deg/ms and 
0.26+0.04 deg/ms for Rat! and Rat5, respectively. For 
repulsion EMF, the slopes of eye movement-current duration 
are 0.44+0.05 deg/ms and 0.1+0.03 deg/ms for Rati and Rat5, 
respectively. 
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Figure 4 Effect of EM coil current intensity on EMF-induced 
eye movements 

A: Representative averaged eye movement responses of Rat2 at two current 
intensities for the two polarities. Error bars indicating standard error are 
within the thickness of the lines. B: Summary of amplitudes of EMF-induced 
eye movements as functions of current intensity from 4 rats with attraction 
EMF (black symbols and lines) and 3 rats with repulsion EMF (grey symbols 
and lines). EM distance from the magnet on the rat eye (10 mm) and EM 
duration (20 ms) remained unchanged throughout the measurements. Error 
bars indicating standard errors are within the thickness of the symbols. 


Effects of the magnet-metal rod distance on the 
EMF-induced eye movement 

Since the EMF is inversely related to square of distance 
between the magnet and the magnetized metal rod, the 
amplitudes of the EMF-induced eye movements were expected 
to be inversely related to square of the distance. Figure 6 
plotted the amplitudes of EMF-induced eye movements as a 
function of the distance for two rats. The two solid lines show 
that the following equations provide excellent fits of the data. 


E = 1424.9/D? (3) 
(R? = 0.99, Rat1) 


E = 786.7/D? (4) 
(R? = 0.99, Rat2) 
Where E is the amplitude of EMF-induced eye movement and 


D is the distance between the magnet and the metal rod 
(Figure 1). 


Main sequence of EMF-induced eye movement 

To assess the dynamics of the EMF-induced eye movements, 
we compared their main sequences to that of saccades in 
awake rats and the eye movements evoked by microstimulation 
of the abducens nucleus in anesthetized rats. Peak velocities 
of the EMF-induced eye movement were plotted against their 
amplitudes for three rats (Figure 7, green/blue/red symbols). 
Data from the three animals fell on a single straight line with 
a slope of 89+1 deg/s/deg (R?=0.98), indicating a well-defined 
main sequence for this condition. However, since EMF had 
a much more rapid onset than the muscle forces generated 
during saccade and abducens nucleus stimulation, the main 
sequence of the EMF-induced eye movements exhibited a 
steeper slope than that of the awake saccade (Figure 7, 
black symbols, slope of 65+10 deg/s/deg, R?=0.62), as well 
as eye movements induced by abducens nucleus stimulation 
in anesthetized rats (Figure 7, grey symbols, slope of 3642 
deg/s/deg, R?-0.9) (P«0.01). 
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Figure 5 Effect of EM current duration on EMF-induced eye 
movements 

A: Representative averaged eye movement responses (Rat1) at three current 
durations for the two polarities. Error bars indicating standard error are within 
the thickness of the lines. B: Amplitudes of EMF-induced eye movements as 
functions of current duration from Rati and Rat5. Other parameters (EM 
distance from the magnet on the rat eye: 10 mm; EM current: 0.55 A) 
remained unchanged throughout the measurements. Error bars indicating 
standard errors are within the symbols. 


Effects of passive stretching on the extraocular muscle 
morphology 

The eyes used in these experiments were stained with Masson 
trichrome stain to evaluate whether repeated FMF-induced eye 


movements resulted in morphological changes. Since the 
EMF-induced eye movements were always in the abduction 
direction for this experiment, any strain was presumed 
to be induced on the medial rectus muscle, which was 
stretched without the tonic eye position motor commands 
being suppressed. The lateral rectus, would be shortened 
under these procedures. Figure 8 shows photomicrographs 
of a section taken through a left globe containing both 
the stretched medial rectus muscle (Figure 8A, B) and 
the compressed lateral rectus muscle (Figure 8C, D). The 
EMF-induced eye movements did not appear to displace 
the medial rectus muscle from its normal location between 
lobes of the Harderian gland (Figure 8A). The medial rectus 
muscle scleral insertion appeared intact. No obvious tearing 
or separation of the muscle insertion were noted. The 
morphology of the neurovascular bundles (not illustrated) 
in the stretched medial rectus muscle was similar to the 
control lateral rectus muscle. Medial rectus and lateral 
rectus muscle structure (boxes in Figure 8A, C) are shown at 
high magnification (Figure 8B, D). The stain revealed similar 
striation patterns for the medial rectus (Figure 8B) and lateral 
rectus (Figure 8D) muscles. No obvious structural changes 
were seen when viewed at higher magnification. Taken 
together, the medial rectus muscle morphology, as determined 
by this qualitative assessment, appeared unchanged by this 
perturbation. We hypothesized that damage to the medial 
rectus muscle structure might manifest as increases in local 
inflammatory response including an increase in inflammatory 
cells or extravasation of blood into the surrounding connective 
tissue. However, no obvious increase in the inflammatory cells 
or accumulation of blood within the muscle or surrounding 
connective tissue was observed in the orbit. Thus, based on 
the qualitative light microscopic analysis, the EMF-induced eye 
movements observed in this study did not appear to produce 
gross changes to the EOM morphology. 


25 


20 


15 


10 


Eye movement amplitude (deg) 





8 10 12 14 


Distance (mm) 


Figure 6 Effect of the distance between the tip of the metal rod 
(C in Figure 1) and the magnet on the rat eye (B in Figure 1) on 
EMF-induced eye movements 

Red dots are for Rat1 and blue dots are for Rat2. EM current was 0.55 A and 
EM duration was 20 ms. Error bars indicate standard errors are within the 
symbols. 
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Figure 7 Main sequence of EMF-induced eye movements (each 
rat is indicated by different colored symbols) 

The grey symbols and the black symbols are for abducens stimulation-induced 
eye movements in anesthetized rats and saccades in awake rats, respectively. 
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Figure 8 Morphologic evidence against EMF-induced damage 
in extraocular muscle fibers 

Photomicrographs of Mason trichrome stained sections through the left 
globe at the level of the medial rectus (stretched) and the lateral rectus 
(compressed) muscles. Low power photomicrographs show the organization 
of the medial rectus (A) and lateral rectus (C) muscles. When the medial 
rectus muscle (box in A) and the lateral rectus muscle (box in C) are viewed 
at higher magnification in B and D respectively, they appear the same. Scale 
bars: 1.0 mm for A and C; 0.1 mm for B and D. 


DISCUSSION 


The feasibility of using EMF to generate eye movement was 
first demonstrated by Bohlen & Chen (2016) in a ping-pong ball 
model. In the present study, we presented further evidence 
in a rat model that the EMF induced well-defined passive 
eye movements over a range of amplitudes (up to 15 deg) 
and peak velocities (up to ~1500 deg/s) without causing 
damages to the extraocular muscles. The EMF approach 
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offers important features that make it an excellent tool to 
study neural processing of EOM proprioception signals in 
awake animal models. First, the EMF generates passive eye 
movement without physically contacting the eye. This feature 
is essential for perturbing fixation or ongoing eye movements 
in awake and behaving subjects. Second, the EMF induces 
eye movement with a very short latency (<5 ms). Since a 
typical saccadic eye movement lasts only 20—50 ms, the rapid 
activation and deactivation of the electromagnet minimizes 
potential biases of subjects’ attention and anticipation. Third, 
by selecting appropriate current parameters (polarity, intensity 
and duration), passive eye movements can be generated 
with desired direction, amplitude and duration. This feature 
allows generation of a wide range of passive movements to 
fully characterize neural encoding of the EOM proprioceptive 
signals. Fourth, by varying magnet and rod positions, it is 
possible to perturb different EOMs and generate passive eye 
movements in different directions. This would overcome the 
limitations of the vibratory method, which primarily activates 
the inferior rectus muscle, as the tendons of other EOMs are 
difficult to reach via this technique (Han & Lennerstrand, 1999; 
Velay et al., 1997). 


Previous studies on proprioceptive signals of the EOMs 
There is a long history of eye muscle proprioception 


investigation (for review see Donaldson, 2000). Since it is 
difficult to access the primary trigeminal afferent neurons under 
chronic conditions, early studies of extraocular proprioceptors 
primarily relied on examining the effects of sectioning the 
ophthalmic branch of the trigeminal nerve (Guthrie et al., 1983; 
Lewis et al., 1998) or surgically altering EOM insertions to 
alter the position of an eye in the orbit (Steinbach & Smith, 
1981). Although important insights have been gained in 
these studies, their approaches were invasive, irreversible, and 
the results may have been confounded by central nervous 
system compensation. Investigators have also tried to directly 
manipulate proprioceptive sensory signals while subjects 
performing visuo-motor tasks, by pulling on the eye via sutures 
or suction, or by tugging or vibrating eye muscles. For example, 
Keller & Robinson (1971) trained monkeys to fixate a visual 
target with one eye while the other eye was passively pulled 
by force applied to a contact lens. Abducens neurons on 
the side of the pulled eye were recorded. They found that 
the abducens neuron activity exhibited no change in response 
to the passively generated eye movements, suggesting that 
there was no stretch reflex in the lateral rectus muscle. While 
these results were widely interpreted as evidence against a 
role of the EOM afferent signals in eye movement control, it is 
important to note the study's limitation, i.e., it only generated 
slow passive eye movements and the eye was not moving. 
It did not rule out the possibility that stretch reflexes may 
be present if the eye was moving at high velocity. Indeed, 
another study using contact lenses to impede the movements 
of an eye during saccades showed changes in the contralateral 
eye (Knox et al., 2000). Furthermore, in a recent study іп 
anesthetized rats and squirrel monkeys, Dancause et al. (2007) 
manually rotated the eye, while measuring the EMG activity 


in the lateral rectus muscle. By producing passive eye 
movements of ~25 deg within about 100 ms (eye velocity of 200 
deg/s), they observed typical stretch reflexes in these muscles. 

In addition to using surgical and mechanical approaches 
to manipulate the proprioceptive EOM signal, pharmacological 
approaches have also been used to investigate the functional 
roles of EOM proprioceptive signals. For example, O’Keefe 
& Berkley (1991) injected dibucaine (a high-potency paralytic 
agent) in one eye of anesthetized cats to block sensory 
nerve transmission in the eye. They found that the 
injection reduced frequency and amplitude of spontaneous 
eye movement in both the treated and non-treated eyes and 
concluded that proprioceptive sensory signals from one eye 
are used to modulate movements of the other eye. This is 
consistent with the fact that the unilateral trigeminal inputs 
are distributed bilaterally by trigemino-trigeminal projections 
(Warren & May, 2013). 

In contrast to the extensive studies on the behavioral effects 
of manipulating EOM proprioceptive signals, few studies have 
examined the neural substrates of the EOM proprioceptive 
signals. In an earlier study, Fuchs & Kornhuber (1969) reported 
evoked potential responses in cat cerebellum while their eyes 
being rotated at saccade velocities. In an important recent 
study, Goldberg and colleagues discovered a representation 
of proprioceptive eye position signals in the somatosensory 
cortex of behaving monkeys (Wang et al., 2007). To elucidate 
the nature of these eye position-related signals, they injected 
lidocaine into the retro-orbital space of an eye in awake 
behaving monkeys to block all sensory inputs from the eye. The 
injection substantially reduced eye position-sensitive neuronal 
activity in the somatosensory cortex. Their subsequent studies 
of gain fields in parietal cortex suggest that while corollary 
discharge is used for guiding the initial saccades in a memory 
guided series, proprioception plays the more important role in 
subsequent saccades in the series (Sun & Goldberg, 2016; 
Xu et al., 2012). 


Future studies 

In comparison with neural processing of corollary discharge 
signals related to oculomotor commands (Sommer & Wurtz, 
2008), neural processing of EOM proprioceptive signals 
has received much less attention. However, understanding 
the EOM feedback signals and their interactions with 
the feedforward efference copy signals has important 
functional and clinical implications. It was the finding of 
EOM proprioceptive signals in the somatosensory cortex 
(Wang et al., 2007) that reignited interest in elucidating the role 
that EOM proprioceptive signals play for gaze control and 


spatial perception at different levels of the oculomotor system. 


The present study was an effort to develop the EMF approach 
for effective manipulation of proprioceptive signals without 
changing ongoing motor commands. The results demonstrate 
the feasibility of using EMF approach in other species, including 
monkeys. Ongoing studies are directed at improving the EMF 
approach in two ways. One is to develop a ferrous metal 


ring that can be implanted under the conjunctiva of the eyes. 


The other one is to further develop the EM coil and control 


of the activation current so that both we can either generate 
an eye movement with desired speeds and amplitudes or hold 
the eye at any desired position for a period of time. This 
approach will offer an important technical advancement to 
facilitate behavioral and neurophysiological studies in awake 
and behaving subjects. 
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ABSTRACT 


This research aimed to provide evidence of a 
relationship between digit ratio and depression 
status in the cynomolgus monkey (Macaca 
fascicularis). In stable cynomolgus monkey social 
groups, we selected 15 depressed monkeys based 
on depressive-like behavioral criteria and 16 normal 
control monkeys. All animals were video recorded for 
two weeks, with the duration and frequency of the 
core depressive behaviors and 58 other behaviors in 
12 behavioral categories then evaluated via 
behavioral analysis. Finger lengths from the right 
and left forelimb hands of both groups were 
measured by X-ray imaging. Finger length and digit 
ratio comparisons between the two groups were 
conducted using — Student's t-test. In terms of 
the duration of each behavior, significant differences 
emerged in “Huddling” and five other behavioral 
categories, including Ingestive, Amicable, Parental, 
Locomotive, and Resting. In addition to the above 
five behavioral categories, we found that depressed 
monkeys spent less time in parental and rubbing 
back and forth behaviors than the control group. 
Furthermore, the 4th fingers were significantly longer 
in the left and right hands in the control group 
relative to the depressed monkeys. The second-to- 
fourth (2D:4D) digit ratio in the left and right forelimb 
hands was significantly lower in the control group 
than that in the depressed group. Our findings 
revealed significant differences in finger lengths and 
digit ratios between depressed monkeys and healthy 
controls, which concords with our view that relatively 
high fetal testosterone exposure may be a protective 
factor against developing depressive symptoms (or 


Science Press 


that low fetal testosterone exposure is a risk factor). 


Keywords: Finger length; Digit ratio; 
depressive disorder; Cynomolgus monkey 


Major 


INTRODUCTION 


Major depressive disorder (MDD) is a debilitating psychiatric 
mood disorder that affects millions of individuals globally 
(Gelenberg, 2010). Our understanding of the biological basis 
of MDD is poor, and current treatments are ineffective in a 
significant proportion of cases. This likely relates to the lack of 
human and non-human primate research models compared 
with the dominant rodent models of depression, which possess 
translational limitations due to limited homologies with humans. 
Therefore, a more homologous primate model of depression is 
needed to advance our understanding of the pathophysiological 
mechanisms underlying depression and to provide a sound 
basis for conducting pre-clinical therapeutic trials. 

Social stress plays a major role in the pathogenesis of 
depression (Krishnan & Nestler, 2008). In human research, 
depressive patients, especially women, are more likely to 
experience depression after prolonged stress (Sherrill et al., 
1997). The diagnosis of depression in humans is based on 
various scales. In line with the DSM-V (Diagnostic and 
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Statistical Manual of Mental Disorders, fifth edition) diagnostic 
criteria, MDD is characterized by five (or more) of the following 
symptoms: depressed mood, loss of interest, weight change, 
sleep disturbance, psychomotor agitation or retardation, loss 
of energy, feelings of worthlessness, difficulty concentrating, 
and recurrent thoughts of death. These symptoms must 
persist for at least two weeks and must include depressed 
mood and/or loss of interest (Gnanavel & Robert, 2013). In 
determination of primate depressive behavior, the most 
reliable method is through behavioral phenotypes. Because of 
the inability to transfer pressure after being the subject of 
aggression, depressed female cynomolgus monkeys (Macaca 
fascicularis) can experience long-term social pressure. Shively 
et al. (2005) indicated that socially-subordinate female 
cynomolgus monkeys, who are likely weak competitors in 
social environments due to long-term attack and suppression, 
can exhibit similar pathogenetic processes as depression in 
humans. Female cynomolgus monkeys have several 
behavioral and physiological characteristics in common with 
human depression. They tend to spend more time alone than 
their dominant counterparts, exhibit greater vigilance, display 
slumped or collapsed body posture, show diminished interest 
in feeding and sex, and subdued communication and 
reciprocal grooming with others for at least two weeks (Shively 
et al., 2005). Huddling, which is defined as a slumped body 
posture with the head at or below the shoulders during the 
awake state (i. e., when the monkey’ s eyes are opened) 
accompanied by a relative lack in responsiveness to 
environmental stimuli, is used as a behavioral indicator of 
depression and is a core posture reflecting depressed mood 
in monkeys (Felger et al., 2007). Based on these well- 
established criteria, we successfully constructed a naturally 
occurring depression model in macaques (Xu et al., 2015). 
Because of the similarities in complex behavioral and 
psychological processes between macaques and humans, the 
development of naturally occurring animal models of disease, 
whether physical or psychological, is a valuable approach for 
translational research between human studies and induced 
primate models, especially for depression (Capitanio, 2017). 
The ratio of the index to ring finger length is a commonly 
used measure. The male ring finger is longer than the index 
finger, whereas the index finger of the female is basically 
equal to that of the ring finger. The second-to-fourth (2D: 4D) 
digit ratio was first proposed by Manning et al. (1998) to 
predict estrogen and sperm number and prenatal sex 
hormones. Under regulation of the Hox gene, hormone levels 
during embryonic development, especially the level of 
androgen, affect finger development (Goodman & Scamble, 
2001). The 2D:4D ratio is present in a human embryo by the 
seventh week of pregnancy (Lutchmaya et al. 2004). 
Additionally, evidence suggests that in both males and 
females, the digit ratio can act as an indicator of the level of 
testosterone the developing fetus was exposed to, making it a 
useful indirect measure of organizational prenatal hormone 
exposure. Some research indicates that the 2D:4D ratio may 
also be an indicator of perinatal androgen action, whereby 
lower digit ratios predict greater androgen sensitivity (Smedley 
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et al., 2014). Thus, given that depression is a strongly sexually 
dimorphic trait, it is reasonable to expect that the 2D:4D ratio 
may be related to depression (Evardone & Alexander, 2009). 
Furthermore, previous investigations have confirmed that the 
2D:4D ratio is associated with human behavior, such as 
aggression, cooperation, left-handedness, and human disease, 
including breast cancer, dyslexia, infertility, myocardial 
infarction, and autism (Ronalds et al., 2002). For instance, 
Bailey & Hurd (2005) reported that the 2D:4D ratio is related 
to individual aggressiveness, confidence, and competitive 
ability. Furthermore, other developmental disorders, emotional 
behaviors, and negative and affective symptoms in schizophrenia 
are reported to be related to the 2D: 4D ratio (Carre et al., 
2015). Previous studies have also shown that individuals with 
a higher 2D:4D ratio are more likely to suffer from depression. 
For instance, Smedley et al. (2014) found that a higher digit 
ratio is correlated with higher depression scores in females, 
but not males. To date, however, results have been inconsistent. 
For example, using a large sample comprised of 298 college 
students (149 males and 149 females), Bailey & Hurd (2005) 
found that more feminine ratios were associated with higher 
depression in men, but found no correlation between the 2D: 
4D ratio and depression in women. Bailey & Hurd (2005) 
produced the unusual finding of no sex differences in depression, 
although digit ratios did differ in the expected direction. As 
such, it remains unclear to what degree depression and the 
2D: 4D ratio, both characterized by marked sex differences, 
are related (Smedley et al., 2014). 

Females are more susceptible to depression in social 
groups and depression is approximately twice more common 
among women than men (Trivers et al., 2006). Therefore, we 
chose female cynomolgus monkeys as our research targets. 
Based on our previous method (Xu et al., 2015), a total of 15 
depressed female monkeys were selected across 52 enclosures, 
with 16 healthy subjects selected as controls from the original 
population (n=6 012). A strict radiographic procedure was 
used to measure index and ring finger lengths in both forelimb 
hands. This investigation aimed to provide evidence of a 
relationship between the digit ratio and depression status in 
the cynomolgus monkey. 


MATERIALS AND METHODS 


Ethics statement 

Behavioral data acquisition was observational under normal 
circumstances and did not involve physical manipulation of 
the subjects or changes to their environment or diet. Animal 
care and housing procedures followed Chinese regulatory 
requirements and the Association for Assessment and 
Accreditation of Laboratory Animal Care International. In brief, 
complete animal husbandry and veterinary care were 
provided daily. Animals were fed a nutritious standardized diet, 
supplemented daily with fresh fruits and vegetables. Animals 
had unrestricted access to potable water and their enclosures 
were cleaned each day. Animals were observed daily by 
trained care-takers. Any observed abnormality, disease, or 


injury was reported to the veterinary staff for diagnosis and 
treatment; this veterinary support was documented in both 
hard copy and electronic formats. In addition, this study was 
performed in strict accordance with the recommendations in 
the Guide for the Care and Use of Laboratory Animals of the 
Institute of Neuroscience of Chongqing Medical University 
(Approval No.: 20100031). Prior to implementation, the 
experimental protocol was approved by the Committee on the 
Ethics of Animal Experiments at Chongqing Medical University 
and was in accordance with state regulations. 


Observation site 

The M. fascicularis Feeding Base of Zhongke Experimental 
Animal Co., Ltd. is in Suzhou, China, at E31°07'03" to E31? 
07'06", N120?19'08" to N120?19' 15". The company imported 
the M. fascicularis subjects from Guangdong Province, China 
and from Vietnam in 1990, from which they established a 
domestication and breeding base for these monkeys. 


Subjects 
We scanned a total population of 6012 adult female 





cynomolgus monkeys across all 52 enclosures. Depressive 
behavior was identified using the operational definition 
according to Shively' s criteria: slumped or collapsed body 
posture (Figure 1), diminished interest in feeding and sex, and 
diminished communication and reciprocal grooming with 
others (Xu et al., 2015). Sixteen healthy adult female M. 
fascicularis subjects (aged 9—13 years) were randomly 
selected from the pool of 6 012 monkeys. A total of 15 
depressed female monkeys (aged 10—12 years) were selected 
from the 52 enclosures based on the above-mentioned 
depression phenotypes lasting for at least two weeks. All 
subjects were reared in socially-stable colonies with negligible 
rates of conflict (Willard & Shively, 2012). Staff veterinarians 
ruled out disease in the subjects. Each colony was housed in 
an indoor free enclosure measuring 8.0 mx3.0 mx3.0 m (Lx 
W x H) with continuous daylight exposure. Every colony was 
composed of two males, 16—22 adult females, and their 
offspring of less than six months of age. To reflect wild 
populations, the male:female ratio was maintained at 
1:(7-10). 


Figure 1 Core depressive behavior (“Huddling”) in a cynomolgus monkey 
A: Normal cynomolgus monkey; B: Depressed cynomolgus monkey exhibiting “Huddling” behavior. 


Behavioral recording methods and scored behavioral items 
are described in our previous work (Xu et al., 2012). The 
duration and frequency of the core depressive behavior 
"Huddling" (Figure 1) and 58 other behaviors in 12 categories 
(Ingestive, Thermoregulatory, Rutting and estrous, Mating, 
Resting, Parental, Amicable, Conflict, Vigilance, Communication, 
Locomotive, Miscellaneous behaviors) were video recorded 
by three well-trained observers blind to the behavioral definition 
using NOLDUS Observer XT software (v10.0, Noldus Information 
Technology, Leesburg, PA, USA) during two consecutive 
weeks with four phases per day (A1 0900-0930 h; A2 0930- 
1000 h; P1 1500-1530 h; P2 1530-1600 h). Behavioral data 
were coded as duration (in seconds) and frequency (in count) 
for each discreet behavioral item per 30 min observational 
phase and presented as means+SD. 


Finger length measurement with radiography 
We measured index and ring finger lengths from the right and 
left forelimb hands in both groups after the last day of 


behavioral data acquisition using a digital radiography unit 
with a flat-panel digital detector (PLX8200, Perlove, Nanjing, 
China) (Kalichman et al., 2017). The digital detector was 
exposed to X-ray at 60 kVp, with an approximate detector-to- 
tube distance of 1 m. Exposure times were no greater than 
0.1 s, resulting in 4.0 mA exposure. Three qualified staff 
performed this process in cooperation: one undertook 
anesthesia and hand-position adjustment, one operated the 
machine to acquire finger length, and one recorded the data 
and was blind to the experiment. Ketamine (10 mL / kg) 
anesthetic was administered intramuscularly (I.M.) in the distal 
hind limb at 0900 h (Nelson & Voracek, 2010). Approximately 
10 min after ketamine injection, the finger length ratio was 
measured at an accuracy of 0.01 mm using X-ray imaging 
(Choi et al., 2011) (Figure 2). 


Statistical methods 
To assess the behavioral differences between depressed 
subjects and healthy controls, Student' s t-test was performed 
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Figure 2 Forelimb hand X-ray image 
2D: Index finger; 4D: Ring finger. 


if the data were normally distributed; otherwise, the Mann- 
Whitney U test was applied. As there were 58 behavioral 
items in total, Bonferroni correction was used to reduce type | 
errors for multiple comparisons. Primate fingers are termed 1 
through 5 from the thumb to the smallest digit, respectively. In 
this investigation, the 2D:4D digit ratio was used (Kyriakidis et 
al., 2010). To minimize error, all finger length data were 
presented as теапѕ+ 50. Finger length апа digit ratio 
comparisons between the two groups were conducted using 
Student’ s t-test. P values of less than 0.05 were deemed 
significant for all data. Data management and statistical 
analysis were performed using SPSS 21.0. 


RESULTS 


Differences in behavioral phenotype between depressed 
and control monkeys 

In terms of the duration of behavior, depressed monkeys 
spent more time "Huddling" (P«0.001), in which the target 
displayed a fetal-like, self-enclosed posture with a lowered 
head during the awake state (i.e., when the monkey' s eyes 
were opened). The increased duration of “Huddling” indicated 
that the depressed subjects recursively stayed in the 
depressed mood condition. For Ingestive behavior, the control 
subjects preferred to "feed while sitting" (P«0.05), whereas 
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the depressed subjects displayed more vigilance and 
preferred “feeding while hanging” to avoid potential threats 
and attacks. They also spent more time “licking food residue 
off the cage” (P<0.05) due to the pressure of other monkeys 
during normal food intake. In accordance with the depression 
criteria, depressed monkeys also spent less time drinking. The 
impact of pressure on both groups was further reflected in 
resting behavior, with the depressed group more reluctant to 
“hang on iron chain rest” (P<0.001) compared to the control 
group, but more willing to choose a remote area for resting, 
namely, “hanging on skylight rest” (P<0.001). Regarding 
locomotive behavior, depressed individuals exhibited less 
vitality in walking and standing and spent relatively less time 
performing “quadrupedal walking on floor" (P«0.001), “walking 
on iron chain" (P«0.001), and "standing" (P«0.05). Indicative 
of a friendly relationship among others, the control group 
received more amicable grooming (i.e., “being groomed” (P< 
0.001) and groomed others more often (i. e., “mutual 
grooming” (P<0.001)) than exhibited by the depressed group, 
suggesting reduced interaction in depression. Furthermore, 
there was a reduction in the duration of “nursing” (P<0.05) 
parental behavior. For mating behavior, which is a sign of 
sexual interest, the duration of “copulation” (P<0.05) was 
higher in the control monkeys in comparison to the depressed 
females, with a significant difference in the frequency of 
mating behavior also observed (Table 2). In the end, except 
the frequency of behaviors matched with the duration of 
behaviors, a lower frequency of miscellaneous behavior (i. e., 
"rub hand back and forth" (P«0.001)) existed in depressed 
individuals, suggesting that depressed monkeys may be less 
imposing and lacking in confidence. 


Finger length data 

Finger length was measured in the 16 control and 15 
depressed animals (Table 3). For both the left and right hand, 
the ring finger was significantly longer in the control group 
than that in the depressed group. 


Digit ratio comparison 

The digit ratio was significantly lower in the control group than 
in the depressed monkeys (Table 4), including the 2D:4D ratio 
in the right and left forelimb hands of depressed and control 
subjects. 


DISCUSSION 


With respect to primate finger length investigations, the 2D:4D 
ratio is strongly related to social behavior and physical 
aggression. Nelson & Shultz (2010) reported that a low 2D:4D 
ratio is associated with more competitive social systems, 
which is in accordance with our previous observation (Zhou et 
al., 2014) that depressed monkeys face greater competition 
for social  resources-including feeding ^ opportunities, 
comfortable resting places, and mating opportunities—and 
display significant deficits in social interactions. Thus, the 
cynomolgus monkey population is a suitable choice to study 
the relationship between the 2D:4D ratio and depression. 


Table 1 
Behavior 
Behavioral category 
Core behavior 


Ingestive behavior 


Amicable behavior 


Parental behavior 


Locomotive behavior 


Resting behavior 


Data are means+SD. Bonferroni correction was used to reduce type | errors for multiple comparisons. Mann-Whitney U test, 


*** P«0.001. 


Behavior 
Huddling*** 


*** 


Drinking 


Feeding while hanging* 


Lick food residue off cage* 


*** 


Feeding while sitting 
Mutual grooming*** 
Being groomed*** 


Nursing* 


Quadrupedal walking on floor 


Walking along iron chain 


Standing* 


Hanging on iron chain rest*** 


kkk 


Hanging on skylight rest*** 


kkk 


Duration of behaviors observed in depressed and control subjects 


Duration 
Depressed group 
365.44+514.14 
2.51+13.09 
6.16+33.19 
1.29+13.78 
29.73+89.06 
51.96+153.61 
51.75+120.17 
7.33+36.69 
62.82+58.99 
0.13+1.26 
9.75+25.35 
5.85+71.70 
23.18+143.93 


Table 2 Frequency of behaviors observed in depressed and control subjects 


Behavior 


Frequency 


Control subjects 
159.00+355.08 
5.92+21.59 
3.56+24.14 
0.44+3.24 
63.29+171.94 
114.13+226.53 
97.94+198.63 
11.75+57.42 
82.03+75.68 
0.85+5.62 
13.94+29.95 
45.35+176.59 
8.88+81.59 


P value 
P<0.001 
P<0.001 
P=0.0004 
P=0.0007 
P<0.001 
P<0.001 
P<0.001 
P=0.0006 
P<0.001 
P<0.001 
P=0.0006 
P<0.001 
P<0.001 


Adjusted P value 
P<0.001 
P<0.001 
P=0.0243 
P=0.0419 
P<0.001 
P<0.001 
P<0.001 
P=0.0334 
P<0.001 
P<0.001 
P=0.0378 
P<0.001 
P<0.001 


*: P«0.05; ** P«0.01; 





Behavioral category 
Core behavior 


Ingestive behavior 


Amicable behavior 
Mating behavior 
Parental behavior 


Resting behavior 


Locomotive behavior 


Miscellaneous behavior 


Behavior 
Huddling*** 


Feeding while hanging*** 


Feeding while sitting** 


de 


Drinking 
Mutual grooming*** 
Copulation* 


Nursing* 


Hanging on iron chain rest*** 


Hanging on skylight rest*** 
Walking along iron chain 
Rub hand back and forth*** 


kkk 


Depressed group 
1.53+2.17 
0.11+0.63 
0.85+2.15 
0.20+0.86 
0.74+1.83 
0.08+0.35 
0.28+1.01 
0.02+0.22 
0.29+1.63 
0.03+0.21 
0.02+0.15 


Control subjects 
0.51+1.25 
0.05+0.28 
1.11+2.47 
0.51+1.59 
1.33+2.40 
0.12+0.53 
0.40+1.44 
0.23+0.80 
0.10+0.62 
0.10+0.53 
0.06+0.43 


P value 
P<0.001 
P<0.001 
P<0.001 
P<0.001 
P<0.001 
P=0.0007 
P=0.0002 
P<0.001 
P<0.001 
P<0.001 
P<0.001 


Adjusted P value 
P<0.001 
P<0.001 
P<0.01 
P<0.001 
P<0.001 
P=0.0418 
P=0.0129 
P<0.001 
P<0.001 
P<0.001 
P<0.001 


Data are means+SD. Bonferroni correction was used to reduce type | errors for multiple comparisons. Mann-Whitney U test, *: P<0.05; **: P<0.01; 


m P«0.001. 


Table 3 Digit length in right and left forelimbs of depressed and 


control subjects (cm) 


depressed and control subjects 


Table 4 20:40 ratio іп the right and left forelimb hands of 

















Control Depressed Гей Digit Control Depressed 
ш r: i К P value 
Left/Right Digit: (n716) (n=15) Right Ratio n Mean SD n Mean SD 
Mean SD Mean SD Left 20:40 16 0.79 0.04 15 083 0.03 0.002** 
2 29.76 246 29.28 2.62 Right 2D:4D 16 0.81 0.03 15 0.84 0.05 0.0372* 
Left 4* 37.75 2.60 35.45 3.17 Data are means+SD. Student’ s t-test, * P<0.05; **: P<0.01; ***: P< 
0.001. 
2 30.25 2.67 29.04 2.78 
Right я 
* 37.14 2.64 35.09 3.24 In the present study, we employed a reliable naturally 


t: Fingers are numbered 1 to 5 from thumb to smallest digit. Data are 


теапѕ+50. Student's t-test, *: P«0.05. 


occurring primate depression model to identify depressed 
cynomolgus monkeys. Our data disclosed that the ring fingers 


in both the left and right forelimb hands were longer in healthy 
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female monkeys than in depressed female monkeys. In terms 
of the finger length data, we found that the digit ratios were 
significantly higher in depressed monkeys, including the 2D: 
AD ratios in the left and right forelimb hands. Depression can 
occur due to long-term social pressures, especially in females 
(Sherrill et al., 1997). In a competitive environment, the dorsal 
anterior cingulate cortex (dACC) region in the brain controls 
support, emotion regulation, conflict monitoring, and behavioral 
inhibition (Dedovic et al., 2016). Gorka et al. (2015) revealed 
a significant positive correlation between the 2D:4D ratio and 
gray matter volume of the dACC in women but not in men. 
Interestingly, maturation of the dACC influences the 
development of MDD (Ho et al., 2017). The critical hippocampal 
brain area, which is strongly associated with the pathogenesis 
of depression, is also related to the 2D:4D ratio (Kallai et al., 
2005). These studies provide a possible intrinsic link in the 
brain tissue between 2D:4D ratio and depression. 

Previous investigations have shown that the digit ratio 
persists in a stable range during embryogenesis and increases 
in accordance with personal growth (Goodman & Scamble, 
2001; Herault et al., 1997). Interestingly, Williams et al. (2000) 
used the 2D:4D ratio to reflect the degree of prenatal 
androgen exposure in humans. However, a growing body of 
evidence indicates that the 2D:4D ratio is unrelated to adult 
sex hormone (e.g., estrogen and androgen) concentrations 
(Muller et al., 2011). The 2D:4D ratio appears to be relatively 
stable, although it does increase somewhat throughout 
childhood (Trivers et al., 2006). Thus, there is a general 
consensus that the 2D:4D ratio is a relatively stable biomarker 
for the balance between fetal testosterone (FT) and fetal 
estrogen (FE), with low FT and high FE linked to high 2D:4D 
(Manning et al., 2014). Based on these findings, we hypothesize 
that the 2D and 4D finger length ratios are primarily 
determined by prenatal sex hormone exposure, and that the 
effects of this prenatal hormone on the 2D:4D ratio are not 
presented as estrogen or androgen concentrations differences 
in a later period. However, the prenatal sex hormon affect the 
subject's neural development and biochemistry (Honekopp et 
al., 2007). Those with longer finger lengths tend to possess 
poor aggression tendencies and emotion regulation function. 
This is also consistent with our findings. In a competitive 
environment, high 2D:4D female individuals exhibited a high 
correlation with depression. Thus, our findings provide a novel 
way in which to select depressed monkeys according to 
comparison of the 2D:4D ratio. Future work should examine 
the relationship between the 2D:4D ratio and the severity of 
depression in larger samples that report a wider range of 
depression symptoms. Measurement of the 2D:4D ratio may 
provide a predictive tool for the diagnosis of depression and 
strong support for indications of depression risk to proceed 
early intervention. 


CONCLUSIONS 


Most previous primate digit ratio studies have been examined 
in regard to social behaviors and rank. However, few have 
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investigated the relationship between digit ratios and 
depression in primates. This is the first study to reveal 
significant differences in finger lengths and digit ratios between 
depressed monkeys and healthy controls. We discovered that 
depressed monkeys presented with shorter 4th fingers and a 
higher 2D:4D ratio in both forelimbs. These metrics show 
promise as gross biological indicators to facilitate screening 
for depressed monkeys in large population-based studies. 
However, whether this conclusion can be applied to screen for 
human depression requires further investigation. 
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ZOOLOGICAL RESEARCH 


High egg rejection rate in a Chinese population of grey- 
backed thrush (Turdus hortulorum) 


DEAR EDITOR, 


Several previous studies have indicated that nest sanitation 
behavior is a general adaptation in altricial birds, with egg 
recognition capacity evolving as a specific response to 
interspecific brood parasitism (IBP). However, a recent study 
suggested an alternative hypothesis, concluding that conspecific 
brood parasitism (CBP) selects for egg rejection in thrushes, 
with IBP as a by-product. In the present study, we used a 
spectrophotometer to quantify egg coloration and egg mimicry 
and performed artificial parasitism experiments in the grey- 
backed thrush (Turdus hortulorum). We showed that individuals 
of this species rejected 100% of 12 foreign eggs, without IBP 
or CBP detected. In a review of previous studies, we also 
discuss possible explanations for the high egg rejection rate in 
the grey-backed thrush and suggest areas for future study. 
Altricial birds have evolved advanced reproductive behavior 
to increase the fitness of their offspring by building elaborate 
structures (i.e., nests), in which they lay eggs and rear their 
nestlings (Hansell, 2000). Bird nests not only provide a 
suitable place for nestling development, but also act as a 
concealed location for safety from predators. Furthermore, 
bird parents have evolved nest sanitation behavior to clean 
foreign objects from their nests, including feces, eggshells, 
branches, and leaves, because they induce predation, facilitate 
microorganism growth, damage eggs, or hurt nestlings during 
brooding (Guigueno & Sealy, 2012). Therefore, nest sanitation 
has evolved as a general behavior in altricial birds for 
distinguishing between egg-shaped and non-egg-shaped 
objects. This recognition capacity has further improved in 
some species to facilitate detection of differences within eggs 
(i.e., egg recognition) as a specific adaptation to avian brood 
parasitism, where other birds lay parasitic eggs in nests that 
are not their own, thereby reducing reproductive output of the 
hosts (Davies, 2011; Yang et al, 2015a) Avian brood 
parasitism can be classified as either interspecific brood (IBP) 
or conspecific brood parasitism (CBP). Numerous empirical 
studies have shown that IBP selects for the capacity of hosts 
to recognize and reject foreign eggs (Davies, 2000; Liang et 
al., 2016; Moksnes et al., 1991; Rothstein & Robinson, 1998; 
Yang et al, 2010). Alternatively, the "collateral damage" 
hypothesis states that CBP is responsible for egg rejection in 
birds, with rejection due to IBP constituting a by-product of 
host adaptations against CBP (Jackson, 1998). However, this 
hypothesis failed to explain egg recognition by hosts because 
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it was tested and supported in a single non-passerine bird 
species (Lyon & Eadie, 2004, 2008). Recently, new research 
re-examined this hypothesis and drew supportive conclusions 
by testing it in two passerine species of thrush, that is, the 
song thrush (Turdus philomelos) and European blackbird 
(Turdus merula) (Samas et al, 2014a) However, Soler 
(2014a) stated that, to date, there is no evidence of CBP 
causing egg rejection in thrushes per se, though Samas et al. 
(2014b) subsequently supported their conclusion with 
empirical evidence. Similarly, Ruiz-Raya et al. (2016) 
investigated recognition of conspecific or heterospecific eggs 
in European blackbirds by manipulating the risk of CBP and 
IBP, respectively. They found that blackbirds presented low 
recognition of conspecific eggs even under high risk of CBP, 
and thus their results supported the IBP hypothesis that egg 
recognition has evolved and is maintained in blackbirds as a 
response to previous cuckoo parasitism. 

Here we performed an empirical study to test egg recognition 
capacity in the grey-backed thrush. The main purpose of this 
study was to provide initial information on egg recognition in 
this species, which may facilitate further study. According to 
previous studies on Turdus spp., we predicted that the grey- 
backed thrush would not show egg recognition capacity 
because no CBP or IBP has been found in this population. An 
alternate prediction was also considered, that the grey-backed 
thrush may also display egg recognition due to previous IBP 
by cuckoos, which still affects host behavior. 

This study was performed in Fusong County, Jilin Province, 
China (N42? 19' 382", E127? 15' 107"), an area of secondary 
forest fragmented by corn (Zea mays) crop farmland and 
scattered plantations (dominated by larch Larix spp.), from 
May to June 2013. This region is in the temperate zone at an 
elevation of 481 m, with a continental monsoon climate 
characterized by cold and snowy winters with an average 
annual temperature of 4 ?C. The grey-backed thrush belongs 
to the Turdidae family and is mainly distributed in East Asia 
(MacKinnon & Phillipps, 1999), where it chooses nest sites 
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with short ground cover and a high density of small trees and 
shrubs (Zhou et al., 2011). In our study site, open-cup nests 
were built in trees (Figure 1A), and pale green eggs with 
reddish markings were laid in these nests. (Figure 1B), with an 
average clutch size of 4.42 eggs+0.51 (range 4—5 eggs, n=12). 





Figure 1 Nest site, nest, incubating female, and eggs of the grey- 
backed thrush (Photos by Long-Wu Wang) 


All experiments complied with the current laws of China. 
Experimental procedures were in accordance with the Animal 
Research Ethics Committee of Hainan Provincial Education 
Centre for Ecology and Environment, Hainan Normal University 
(permit No. HNECEE-2012-004). 

The appearances of the thrush and model eggs were 
quantified with a spectrophotometer (Avantes-2048, Avantes, 
Apeldoorn, The Netherlands). Six reflectance spectra were 
measured in each egg and averaged to represent the color of 
the egg. Model eggs were immaculate blue, acommon coloration 
of cuckoo eggs in China (Yang et al., 2010, 2012). Thrush 
eggs are pale green with reddish markings and thus their egg 
ground color and markings were measured, respectively. 
Subsequently, egg spectra were loaded into AvaSoft 7.3, at 
which the wavelength range of spectra varied from 300 to 
700 nm. The spectral range of 300—400 nm and 400-700 nm 
refers to ultraviolet (UV) light and visible light (VIS), 
respectively (for details, see Yang et al., 2010, 2013). 

Nests of the grey-backed thrush were found by searching 


potential nest sites and monitoring the activity of reproductive 
adults. Nests were then randomly sorted into two groups: (1) 
manipulated group into which one immaculate blue model egg 
was introduced just after clutch completion (n=12) (Figure 2); 
and (2) control group in which nests were visited by the same 
procedure to control for human disturbance, but no manipulation 
was made (n=10). Manipulation was performed in these 
circumstances without observing hosts to avoid potential 
effects of host observations on recognition (Hanley et al., 
2015). Observed nests were monitored for 6 d after manipulation 
and the responses of thrushes to foreign eggs were classified 
as rejection, if foreign eggs were ejected, pecked, or deserted, 
or accepted if foreign eggs were intact and incubated (Yang et 
al., 2010, 2014b). Model eggs were made by polymer clay 
and their sizes were standardized to 25 ттх 19.5 mm, similar 
to thrush eggs (25.07+0.42 ттх19.74+0.71 mm, n=10). 





Figure 2 Experimental nest of the grey-backed thrush with a blue 
model egg (Photo by Long-Wu Wang) 


The grey-backed thrushes laid pale green eggs with dense 
reddish markings (Figure 1B). Egg reflectance analysis 
illustrated that egg ground color was consistent with human 
eye assessment, with a reflectance peak in the range of green 
light (475—550 nm) (Figure 3). Similarly, the reddish markings 
had a reflectance peak in the range 550—620 to 620-700 nm, 
which represents yellow and red light, respectively (Figure 3). 
The blue model egg had a reflectance peak in the range of 
blue light (400—475 nm) (Figure 3). The reflectance contrast 
between trough and peak reflects chroma (or color saturation). 
Therefore, blue model eggs were more saturated in color than 
the thrush egg. In brief, the model egg was very different from 
thrush eggs according to vision based on human eyes and 
spectral reflectance. Experimental parasitism indicated that 
the grey-backed thrush rejected all non-mimetic model eggs, 
with a rejection rate of 10096 (n=12). All rejections occurred 
within a day (i.e., 24 h) and were all performed by ejection 
without any recognition error. No rejection was found in the 
control group (n=10). 

This Chinese population of grey-backed thrush possessed 
high recognition capacity, rapidly rejecting 100% of non- 
mimetic foreign eggs. This contradicted our expectation that 
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Figure 3 Egg reflectance of the grey-backed thrush and model eggs 
Reflectance of thrush eggs was averaged from three eggs from the 
nest in Figure 2. 


grey-backed thrush should have no or low-level egg 
recognition capacity because no CBP or IBP was detected in 
our study population. Thus, this egg rejection ability of the 
grey-backed thrush needs further investigation. 

A recent study on artificial parasitism in song thrushes and 
European blackbirds found unusually high rejection rates of 
CBP (up to 60%) and positive co-variance with conspecific 
population densities without risk of IBP (Samas et al., 2014). 
Because IBP rejection rates did not covary geographically with 
IBP risk (Grim & Stokke, 2016) and thus contradicted the IBP 
hypothesis that egg recognition has evolved as a response to 
IBP, Samas et al. (2014) concluded that egg recognition in 
Turdus spp. has evolved as a response to CBP, not to IBP. 
However, if egg rejection abilities can be maintained in the 
absence of IBP, high egg rejection rates can still be retained 
without geographic co-variation. Therefore, Soler (2014a) 
criticized the conclusion of Samas et al. (2014) and argued 
they made an invalid conclusion due to an out-of-date 
theoretical background and a biased selection of references. 
However, Samas et al. (2014b) argued that a theory is never 
out of date and addressed the theoretical objections by 
empirical evidence. Recently, Ruiz-Raya et al. (2016) further 
tested egg recognition in blackbirds by manipulating the risk of 
CBP and IBP and concluded that selection from IBP likely 
accounts for egg recognition in blackbirds. 

In our study population, grey-backed thrushes displayed 
high recognition capacity of foreign eggs. It is generally 
accepted that IBP rather than CBP contributes to egg 
recognition in hosts. Firstly, egg recognition capacity is much 
more unlikely to evolve in response to CBP than IBP because 
IBP gives rise to dramatic fitness costs, which are much lower 
than those from CBP (Petrie & Møller, 1991). Parasites and 
hosts with CBP are conspecific and share the same gene pool 
and thus constitute much weaker selection than IBP (Ruiz- 
Raya et al., 2016; Soler, 2014b). Furthermore, conspecific egg 
phenotypes in CBP are too similar to initiate evolution of egg 
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recognition (Soler et al., 2011). However, intraspecific variation 
in egg coloration is high for some host species, and CBP may 
account for the evolution of egg recognition (Cassey et al., 
2008a, 2008b; Hanley et al., 2017; Samas et al., 2014a). 
Secondly, in order to select for conspecific egg rejection, the 
level of CBP must be high. However, in our study population, 
no CBP was detected. Similarly, Samas et al. (2014) described 
that the rates of CBP were only 0%-2.2% and 0%-3.1% for 
the song thrush and blackbird, respectively. However, current 
parasitism rates may not represent actual selection pressure 
from IBP and CBP without considering other factors, such as 
recognition error and rejection cost. Furthermore, because 
egg rejection capacity has evolved in response to IBP, and it 
can be maintained in the absence of IBP, this may also occur 
in response to CBP but in its absence. Many currently 
unparasitized potential host species exhibit a rejection rate of 
nearly 100% (Lahti, 2006; Peer & Sealy, 2004; Yang et al., 
2014a, 2015b). For example, blackbirds were introduced in 
the nineteenth century to New Zealand, where a high rejection 
rate of non-mimetic eggs has been reported (62%, Samas et 
al., 2014; 83.9%, Hale & Briskie, 2007), similar to the rejection 
rate of 90% in Europe (Grim et al., 2014; Martín-Vivaldi et al., 
2012; Moskát et al., 2003). Moreover, a recent review 
concluded that it is not correct to formulate predictions that 
assume that rejection behavior of hosts must be lost in the 
absence of obligate brood parasites (Soler, 2014b). Finally, 
aggressive behavior towards adult cuckoos and reluctance to 
feed cuckoo chicks has been empirically shown in thrushes 
(Grim et al., 2011), providing evidence for contact with IBP in 
the past, resulting in successful resistance against interspecific 
brood parasitism (Ruiz-Raya et al., 2016). However, Samas et 
al. (2014b) argued that aggression in blackbirds did not 
specifically evolve in response to IBP because they are 
aggressive not only to cuckoo dummies but also to any 
intruders near their nests, including harmless pigeons (Columba 
livia). Thus, switching to new types of food is an unlikely 
defense against brood parasite chicks because such changes 
would not prevent most costs from IBP (Grim et al., 2011; 
Samas et al., 2014b). However, Ruiz-Raya et al. (2016) found 
that blackbirds were able to recognize and eject heterospecific 
eggs at high rates, whereas most conspecifics eggs were not 
recognized. Moreover, ejection rates of conspecific eggs did 
not exceed 13%, even in the presence of a high risk of CBP, 
whereas ejection rates of experimental eggs simulating IBP 
were much higher (80%—100%). Female blackbirds were also 
found to be more aggressive towards cuckoos than towards 
blackbird dummies (Ruiz-Raya et al., 2016). Additionally, Ruiz- 
Raya et al. (2016) estimated that the level of CBP necessary 
to select for evolution of host response against conspecific 
eggs would range from 55% to 65%. Because the grey- 
backed thrush has retained a high level of egg recognition, the 
rejection costs, which would contradict such maintenance, 
should be negligible. According to our results, no rejection 
cost was detected. Like blackbirds, the grey-backed thrush 
also has a large bill to grasp foreign eggs for rejection. 
Therefore, rejection costs should not prevent grey-backed 
thrushes from retaining egg rejection capacity. 


In summary, previous studies on thrush hosts have provided 
inconsistent conclusions. This situation is complicated by the 
explanation for one species not being suitable for another. Our 
study provides preliminary information, and thus cannot offer 
sufficient evidence to support either the IBP or CBP 
hypothesis. However, considering that previous studies have 
provided strong evidence that hosts affected by IBP can retain 
egg recognition capacities after long-term escape from cuckoo 
parasitism (Lahti, 2006), egg rejection capacity in grey-backed 
thrushes may be maintained because parasitic cuckoos have 
exploited this potential host in the past. Although theoretically 
hosts may also retain egg recognition due to previous CBP, 
like IBP, empirical studies are currently insufficient to clarify 
this assumption. Furthermore, recent studies have revealed 
that egg accepters can become rejecters after stimulation, in 
all cases switching from acceptance to rejection, implying that 
historical cases of IBP or CBP are considerably underestimated 
in currently non-parasitized potential host species (Molina- 
Morales et al., 2014; Yang et al., 2015b). Further studies 
referring to egg color variation and recognition with different 
degrees of egg mimicry are necessary in the grey-backed thrush 
and even other species of Turdus. Such studies will help us 
better understand the origin of egg recognition in thrushes. 
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Threshold for maximal electroshock seizures (MEST) 
at three developmental stages in young mice 


DEAR EDITOR, 


Early brain development after birth is extremely dynamic, 
suggesting that potential functional changes occur during this 
period. In this study, the maximal electroshock seizure 
threshold (MEST) was used to explore the electrophysiological 
variation among three developmental stages in young mice 
(no more than 5 weeks old). The induced electroshock seizure 
(ES) behavior of early postnatal mice (1—2-weeks old) differed 
from that during weaning (3 weeks old) and early puberty (4— 
5-weeks old) Thus, we further explored their respective 
characteristic responses to the ES parameters. When the 
stimulation current (SC) was limited to 4.0 mA, only the 1—2- 
week-old mice were induced to exhibit ES behavior at 
voltages of 30 V and 40 V, indicating they were more sensitive 
to maximal electroshock seizure (MES) (response to lower 
voltage). Surprisingly, however, they showed substantially lower 
mortality than the older groups under higher voltage 
conditions (60, 100, 160, and 200 V), suggesting better 
tolerance to the SC. We also found that when the current limit 
decreased to 3.5 mA, the 4—5-week-olds mice exhibited stable 
ES behavior with low mortality, while for 3-week-olds mice, the 
SC limit required to be reduced to 1.5 mA. In conclusion, our 
findings showed that neural sensitivity to MES was significantly 
different in young mice before puberty. Thus, greater attention 
should be given to distinguishing the developmental period of 
mice, especially in electrophysiological examination. 

While the macroscopic layout of the brain is nearly complete 
by the end of pregnancy, it develops continuously at a high 
speed until prepuberty (Stiles & Jernigan, 2010). Recent 
magnetic resonance imaging studies have depicted structural 
change processes by brain volume (Gogtay et al. 2004; 
Knickmeyer et al., 2008) and fiber connection (Chen et al., 
2016; Li et al., 2015) growth curves, suggesting that early 
brain development after birth is extremely dynamic (Li et al., 
2013). In addition, many developmental mental disorders likely 
originate from developmental problems in preadolescents 
(Chen et al, 2017; Li et al., 2016). For example, it is 
estimated that 10.5 million children under 15 have active 
pediatric epilepsy, which is more than 10 times greater than 
that found in adults (Keezer et al., 2016). 

Maximal electroshock seizure (MES) is an experimental 
paradigm that induces synchronous neural discharges in the 
brain through artificial current input (Kamei et al., 1978), and 
is used to induce acute epileptic behaviors (Fischer & Muller, 


Science Press 


1988). However, there are limited reports on the application of 
MES in young mice to mimic childhood epilepsy. In this study, 
a MES threshold paradigm was applied to investigate 
electrophysiological variations in young mice less than five 
weeks old. Group information was summarized in Table 1. 
Behavioral expression in the mice included electroshock 
seizure (ES), death, and no response, which reflected their 
brain network states. Young mice less than 5 weeks old were 
divided into three groups: 1—2-week-, 3-week-, and 4—5-week- 
old groups, based on their physical features and activities. 
These groups corresponded to three critical developmental 
stages: i.e., early postnatal, weaning, and early puberty (Miao, 
1997; Xu et al., 2011). 

Usually, an electrical stimulus of the MES paradigm is 
delivered to adult mice (more than 6—week old) and is about 3— 
10 times higher than the individual electrical seizure threshold 
of the animal (Kamei et al., 1978; Murakami et al., 2007). 
The typical MES behaviors are: hind-limb extension, fall, 
and back rigidity, followed by foaming at the mouth and 
urinary incontinence (Ferraro et al., 1998, 2002)(Figure 1A). 
We found that seizure in 3-week-old and 4—5-week-old 
mice induced typical MES behavior as Figure 1A shown. 
Under the same conditions (i.e., SC limited to 3.5 mA and 
voltage of 80 V), the induced seizures of early postnatal 
mice (1—2 weeks old) were different; although urinary 
incontinence and mouth foaming were also observed, the 
seizures of 1—2-week-old mice did not include hind-limb 
extension and fall, with most limbs bending and convulsing 
(Figure 1B). 

These behavioral differences may be due to the sensitivity 
differences in nerves. Therefore, the mice were subjected to 
electrical stimuli with elevated parameters. Surprisingly, when 
the SC was limited to 4.0 mA and the voltage was 160-200 V, 
all 1—2-week-old mice survived the elicited ES behaviors. 
Under the above parameters, all behaviors exhibited by 1—2- 
week-old mice were consistent with typical opisthotonus, 





Received: 30 August 2018; Accepted: 29 January 2019; Online: 28 
March 2019 

Foundation items: This study was supported by the National Natural 
Science Foundation of China (81403191) and Yunnan Provincial Natural 
Science Foundation (2018FB118, 2015FA004, and KKSY201626001) 
DOI: 10.24272/j.issn.2095-8137.2019.038 


Zoological Research 40(3): 231-235, 2019 231 


Table 1 Animals were grouped based on age and current limit in the MEST test 


Current limited to 4.0 mA Current limited to 3.5 mA Current limited to 1.5 mA Total (n) 











1-2 weeks 85 10 10 105 

3 weeks 36 25 52 113 

4-5 weeks 39 70 23 132 

Total (n) 160 105 85 350 
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Figure 1 Electroshock seizure (ES) behavior at different developmental stages in mice 
A: ES behavior of 3-week-old and 4—5-week-old mice. B: Typical ES behavior patterns of 1—2-week-old mice at low voltage (30—100 V, 4.0 mA); C: 


ES behavior patterns of 1—2-week-old mice at high voltage (120—200 V, 4.0 mA). 


which differs from the previous stimulus conditions. The MES under the same conditions (i.e., 4.0 mA and 160—200 V), all 
behaviors in 1—2-week-old mice included bilateral forelimb but one of the 3—week-old and 4—5-week-old mice died after 
clonus, tail stiffness, and bending (Figure 1C). In contrast, rigidity seizures (Figure 2). 
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Figure 2 Maximal electroshock seizure (MES) occurrence, no-response, and death voltage range of 1—2-week-old, 3-week-old, and 4—5- 


week-old mice at 4.0 mA 
When the maximum output current was limited to 4.0 mA, death, MES, and no response were in response to constant voltage electrical stimulation. 


ES: Electroshock seizure. 


Based on the above results, a series of stimuli was given to As shown in Figure 2, when the SC was limited at 4.0 mA 
the three groups of young mice to determine the optimal and the voltages were 30 V and 40 V, ES behavior was 
range of parameters for inducing ES behavior. observed in 1—2-week-old mice (87.5%, n=16) but not in the 
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other groups, indicating that the 1—2-week-old mice were 
more susceptible to electrical stimulation than the older 
groups (Figure 2). This age is equivalent to the human 
breastfeeding stage before 2 years old. Furthermore, this 
result is consistent with the clinical epidemic status of 
epilepsy, a brain disease with abnormal synchronous neural 
discharge of the cerebral cortex, thus indicating that the early 
postnatal brain is more sensitive to electrical stimulations. 

Interestingly, at voltages of 60, 100, 160, and 200 V (Table 
2), 100% of 1—2-week-old mice survived after ES behavior 
occurrence (63/63). In contrast, 95.6% of 3-week-old and 4—5- 
week-old mice died (44/46), demonstrating a substantially 
higher mortality than the youngest group (P«0.000 1). This 
indicated that early postnatal mice better endured electrical 
stimulation than early puberty mice, which has not been 
reported previously. 

Electrical stimuli with reduced current were applied to 
explore the appropriate stimulus range for 3-week-old and 4— 
5-week-old mice. As shown in Figure 3A, when the SC was 
limited to 3.5 mA, 72% (38 in 53) of 4—5-week-old mice 
elicited ES at voltages of 80 V and 100 V. This was a relatively 
safe stimulation range in which to induce ES behavior of 4—5- 
week-old mice. For 3-week-old mice, however, only 50% (4 in 


Maximum current limited at 3.5 mA 


Table 2 Responses of 1 —2-week-old and 3 — 5-week-old mice 
under 4.0 mA and 60-200 V 


Age Number(n) ES Death Chi-square test 
1-2 weeks 63 100% 0 

P«0.000 1 
3—5 weeks 46 4.3% 95.6% 


Statistics were analyzed by Chi-square test. ES: Electroshock seizure. 


8) exhibited successful induction of ES behavior at 60 V. 
Furthermore, when the voltage increased to 80, 100, and 120 V, 
all 3-week-old mice died. When the SC was limited to 1.5 mA, 
3-week-old mice exhibited ES behavior within a broad range 
of voltages (100 – 140 V). This was therefore considered a 
suitable range of the stimulus parameters, although 4—5-week- 
old mice demonstrated no responses under these conditions. 
Comparing the three groups, we determined that the 
behavioral responses of young mice in different periods under 
a MEST paradigm were different due to the different intensities 
of brain network connections. The stimulated parameter 
ranges eliciting MES for different groups of young mice were 
not continuous. Thus, based on our results, it appears that 3 
weeks of age may be a turning point in nervous system 
development in mice. 


Maximum current limited at 1.5 mA 
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Figure 3 Voltage range of 3-week-old and 4—5-week-old mice showing successful maximal electroshock seizure (MES), no-response, and 


death with current limited to 3.5 mA or 1.5 mA 


When the maximum output current was limited to 3.5 mA (left) or 1.5 mA (right), death, MES, and no response of 3-week-old (circle) and 4—5-week- 
old mice (solid circle) were in response to constant voltage electrical stimulation. ES: Electroshock seizure. 


To validate that the behaviors induced by electrical 
stimulation in all young mice were the same as epileptic 
seizures, an antiepileptic positive drug inhibitory experiment 
was conducted. Phenobarbital sodium, which is a commonly 
used treatment for epilepsy, was used here to inhibit the 
epidemiology of MES under an adult median effective dose 
(ED,,) of 2.0 mg/kg. Results showed that it significantly 
inhibited the onset of ES behavior in 1—2-week-old (10/16), 3- 
week-old (7/13), and 4—5-week-old mice (4/8), as shown in 
Table 3. These results suggest that although the behavior and 
ranges of voltage and current that induced ES in young mice 
differed from that observed in adults (Ferraro et al., 2002), 


they were still indicative of epileptic episode and could be 
inhibited by antiepileptic drugs. 

The present study used the MEST paradigm to explore 
neurophysiological differences in three developmental stages 
before puberty in mice. We found that seizure was induced in 
all three groups of young mice (less than 5 weeks of age) 
soon after electrical stimulation. The 1—2-week-old mice 
exhibited both forelimb and tail stiffness with mild convulsions, 
whereas the other two groups displayed hind-limb extension, 
fall, and back rigidity. All young mice foamed at the mouth 
during seizure and experienced urinary incontinence, but 
recovered to normal activity after 2—5 s. In addition to 
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Table 3 Inhibition of phenobarbital sodium on MES behavior 
Total (MES/Anti- 











— 

Group Treatment convulsant) Inhibition (%) 
ES 9 (9/0) 

1-2 weeks 62.5 (P*-0.008 8) 
ES+PB 16(6/10) 

4 ° ES 12(12/0) 53.8 

weeks 

ES+PB 13(6/7) (P=0.005 2) 
ES 16(16/0) 50 

4—5 weeks 
ES+PB 8(4/4) (P=0.006 6) 


ES: Electroshock seizure. PB: Phenobarbital sodium, which was 
administered via intraperitoneal injection at a dose of 2.0 mg/kg. *: 
Statistics were analyzed by Chi-square test. 


differences in seizure behavior, the three age groups 
demonstrated different behavioral outcomes under electrical 
stimulation. For example, 3-week-old mice showed significantly 
higher mortality when the SC was 4.0 and 3.5 mA. Only when 
the current was limited to 1.5 mA could they survive after 
induction of MES, whereas, the mice in the 4-5-week group 
showed no responses to electrical stimulations. The optimal 
MES stimulation conditions in 4—5-week-old mice in this study 
are similar to those of adult rats in previous studies (Ferraro et 
al., 2002). Unexpectedly, we found that not only did 1—2-week- 
old mice experience induced seizures at lower voltage, but 
they also survived at higher voltage stimulation. 

Seizures induced by electroshock are one of the two most 
widely studied animal models of generalized epilepsy, the 
other being pentylenetetrazol (PTZ) administration (Loscher, 
2011). Accumulated evidence implicates structures of the 
brainstem as being involved in both kinds of experimental 
seizures. Stimulation of the midbrain reticular formation induces 
motor seizures in cats, rats, and rabbits (Kreindler et al., 
1958). The substantia nigra also seems to play an important 
role in mediating seizure discharge release. The seizures 
(tonic hindlimb extension behaviors of 3- and 4—5-week-old 
mice in the present study) induced by MES can be prevented 
by lesions of the substantia nigra in rats (Garant & Gale, 
1983). Furthermore, injections of GABA agonist muscimol and 
opiates into the same region can also prevent MES-induced 
seizures (ladarola & Gale, 1982). However, lesions in the 
mesencephalic reticular formation can antagonize the 
production of PTZ-induced convulsions (Jinnai et al., 1969). 
Moreover, seizures induced by PTZ administration rather than 
MES can be protect by bilateral diencephalic lesions. (Mirski & 
Ferrendelli, 1986). It is not yet clear whether the diencephalon 
and substantia nigra are parts of a single complex 
neuroanatomical network mediating experimental seizures or 
whether they belong to two separate independent pathways 
for propagation of different types of seizures. As mice reach 
sexual maturity at 6—7 weeks of age, these young mice 
possibly correspond to different human developmental periods: 
1-2 weeks corresponds to a breastfeeding infant period, 3 
weeks corresponds to the childhood weaning period, and 4—5 
weeks corresponds to early puberty (Miao, 1997). In the new 
postnatal brain, the neural fiber connection is involved in rapid 
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synaptic formation and redundant cutting. For 1—2-week-old 
mice, the connection state ofthe diencephalicand mesencephalic 
structures differs from that of the older group, which may be 
the cause of specific seizure phenotypes and greater resilience 
to electrical stimulation. Previous investigations have found that 
adult mice die from MES due to respiratory arrest (Buchanan 
et al., 2014). The high mortality of the 3-week-old mice may 
also be due to changes in the state of connection in the 
brainstem, resulting in increased susceptibility to inhibition of 
the respiratory center. 

As pediatric epilepsy is one of the most vulnerable diseases 
in children (Keezer et al., 2016), the MEST paradigm was 
used to explore the physiological variation in the prepuberty 
brain. Animal epilepsy models of ES are useful for investigating 
pathophysiological mechanisms and developing or evaluating 
new antiepileptic treatments. However, the development of 
pediatric epilepsy models is still challenging due to dynamic 
development in the immature brain. Yet, MES and intravenous 
pentylenetetrazol tests can be used to demonstrate the 
anticonvulsant properties of anti-epileptic drugs. The former 
acts as an acute seizure model (Castel-Branco et al., 2009; 
Murakami et al., 2007), whereas the latter is a chronic model 
(Loyens et al., 2012). They can help to identify the effects of 
compounds on seizure spread and increase the seizure 
threshold, respectively. In this study, phenobarbital, an 
antiepileptic positive drug, was used to determine if the 
induced behaviors in 1—2-week-old mice were the same as 
epileptic seizures. The MES reactions were compared between 
mice with prior injection of phenobarbital (2.0 mg/kg) and the 
blank solvent control. Results showed that all MES behaviors 
in tested mice at the three developmental stages were 
significantly inhibited at similar intensity with phenobarbital, 
indicating that the MES behaviors of the mice were the same 
as those experienced during epileptic seizures. 

Thus, physiological variation in the prepubertal brain to 
electrical stimulation at different developmental stages was 
found using the MEST paradigm, with the three age groups 
exhibiting different behavioral outcomes. The results from this 
study will improve our knowledge regarding early brain 
development and provide new evidence that neural pathology in 
children differs from that in adults, suggesting that the 
development of the brain from birth to adolescence is extremely 
dynamic. Thus, itis necessary to strictly divide the developmental 
stages of youth to obtain an appropriate animal model. 
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Why humans have large brains with higher cognitive abilities 
is a question long asked by scientists. However, much remains 
unknown, especially the underlying genetic mechanisms. With 
the use of a transgenic monkey model, we showed that 
human-specific sequence changes of a key brain development 
gene (Primary microcephaly 1, MCPH 1) could result in detectable 
molecular and cognitive changes resembling human neoteny, a 
notable characteristic developed during human evolution. This 
study was published in National Science Review (Shi et al., 
2019). 

Neoteny has been hypothesized to help explain why we 
differ from our close relatives, such as the chimpanzee (Pan 
troglodytes) (Liu et al, 2012). Also called juvenilization, 
neoteny is the delay in or slowing down of physiological 
development of a species (Skulachev et al., 2017). Neoteny in 
humans is the retention of juvenile features into adulthood andis 
exaggerated compared to that in non-human primates. The 
neotenous human brain provides an extended window for the 
plasticity of the neural network, and therefore a longer time for 
learning, which may be crucial for the formation of human 
cognition (Bufill et al., 2011). 

In the past two decades, comparative genomic studies have 
identified many candidate genes that carry human-specific 
sequence changes potentially contributing to human brain 
evolution (Bustamante et al., 2005; Clark et al., 2003). However, 
functional dissection of these candidate genes, particularly the 
phenotypic consequences of human-specific mutations, is 
lacking. As a model organism, transgenic mice (Mus musculus) 
have been used in the study of human brain function; for 
example, the study of Foxp2 that carries two human-specific 
amino acid changes (Enard et al., 2009). However, compared 
to the human brain, the mouse brain is dramatically different in 
size, structure, developmental pattern, and function, and 
therefore is not an ideal model to study human brain evolution. 
In contrast, the rhesus monkey (Macaca mulatta), an Old- 
World primate species widely used in biomedical research, is 
a better animal model due to its high sequence similarity with 
humans (>93% for protein coding genes) (Yan et al., 2011; 
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Zhang et al., 2014), yet relatively large phylogenetic distance 
(725 million years of divergence from humans), which alleviates 
certain ethical concerns (Coors et al., 2010). These advantages 
have been demonstrated in the use of transgenic monkey 
models of human diseases (Chen et al., 2017; Liu et al., 2016; 
Luo et al., 2016; Qiu et al., 2019; Yang et al., 2008; Zhang et 
al., 2018). 

MCPH1 is a key gene regulating brain development in humans, 
and its truncate mutations can cause primary microcephaly 
(Jackson et al., 2002). Similarly, MCPH1 knockout in mice and 
monkeys can recapitulate the microcephaly phenotype (Gruber 
et al., 2011; Ke et al., 2016), highlighting its critical role in 
mammalian brain development. In previous research, we 
demonstrated accelerated protein sequence changes in MCPH1 
during primate evolution, especially in the lineage leading to 
the origin of our own species (Wang & Su, 2004). Furthermore, 
in vitro functional analysis showed that human-specific MCPH1 
changes resulted in functional divergence between human 
and non-human primates (Shi et al., 2013). Such lines of 
evidence raise the question as to what phenotypic changes 
would appear following the introduction of human MCPH1 
copies into the rhesus monkey genome. 

Using lentivirus delivery, we successfully generated 11 
transgenic monkeys carrying 2—9 copies of human МСРНТ. 
Brain image analysis showed that the transgenic monkeys 
had similar brain volume and cortex thickness as the controls, 
but their relative brain volume and gray matter percentage 
were higher. Importantly, the transgenic monkeys exhibited 
delayed cortex gray matter development compared to the 
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controls. Furthermore, the transgenic monkeys had more 
immature and fewer mature neurons and glial cells than the 
controls, which, according to the transcriptome data, was 
likely caused by the suppression of many neuron maturation- 
and neuron differentiation-related genes. This observation 
was further illustrated by detailed analysis of developing brain 
laminas at the neurogenesis peak during fetal development. 

Finally, to test whether the observed delay in brain 
development had any impact on behavior and cognition, we 
tested working memory in the transgenic monkeys using 
delayed matching-to-sample tasks. Notably, the transgenic 
monkeys exhibited better working memory and shorter 
reaction time than the controls, suggesting that neotenous 
brain development in the transgenic monkeys was beneficial 
for the formation of cognitive abilities, confirming the proposed 
evolutionary advantage of human neoteny. In our future work, 
we will: (1) examine the molecular pathway explaining how 
human MCPH1 copies delay brain development and neural 
maturation, thus identifying more neoteny-related genes; and 
(2) apply more sophisticated tests to further understand how 
the neotenous brain affects cognitive formation in transgenic 
monkeys. 

This study represents the first attempt to utilize a transgenic 
monkey model to study human brain evolution. It also provides 
the first molecular genetic evidence showing neotenous 
changes during brain development caused by human-specific 
mutations of a single gene. The results highlight the great 
potential of non-human primate models in studying human 
evolution and may pave the way for future studies to explore 
the genetic mechanisms of human-specific traits and elucidate 
the etiology of human brain disorders such as autism and 
Alzheimer’s disease. 
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